Commit d97da3cf authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1292 Rivlin cube ok.

parent d1808580
......@@ -176,7 +176,6 @@
BE01DAD71E854DE300F3EAF7 /* InvariantComputation.hpp in Headers */ = {isa = PBXBuildFile; fileRef = 135F218A1E7ADBC000B9E633 /* InvariantComputation.hpp */; };
BE01DAD81E854DE300F3EAF7 /* InvariantComputation.hxx in Headers */ = {isa = PBXBuildFile; fileRef = 135F218B1E7ADBC000B9E633 /* InvariantComputation.hxx */; };
BE01DADA1E854DE300F3EAF7 /* InvariantHolder.hpp in Headers */ = {isa = PBXBuildFile; fileRef = 135F218D1E7ADBC000B9E633 /* InvariantHolder.hpp */; };
BE01DADB1E854DE300F3EAF7 /* InvariantHolder.hxx in Headers */ = {isa = PBXBuildFile; fileRef = 135F218E1E7ADBC000B9E633 /* InvariantHolder.hxx */; };
BE01DADC1E854DEC00F3EAF7 /* ElementaryDataStorage.hpp in Headers */ = {isa = PBXBuildFile; fileRef = 135F21901E7ADBC000B9E633 /* ElementaryDataStorage.hpp */; };
BE01DADD1E854DEC00F3EAF7 /* ElementaryDataStorage.hxx in Headers */ = {isa = PBXBuildFile; fileRef = 135F21911E7ADBC000B9E633 /* ElementaryDataStorage.hxx */; };
BE01DAE01E854DEC00F3EAF7 /* NumberingSubsetSubMatrix.hpp in Headers */ = {isa = PBXBuildFile; fileRef = 135F21981E7ADBC000B9E633 /* NumberingSubsetSubMatrix.hpp */; };
......@@ -1474,6 +1473,17 @@
BEE45C062357804F00E0AB7C /* Stokes.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE8C36641DB50FB500919468 /* Stokes.cpp */; };
BEE45C072357804F00E0AB7C /* Stokes.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE8C36651DB50FB500919468 /* Stokes.hpp */; };
BEE45C082357804F00E0AB7C /* Stokes.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BE8C36661DB50FB500919468 /* Stokes.hxx */; };
BEE45C092357BA8100E0AB7C /* InvariantHolder.hxx in Headers */ = {isa = PBXBuildFile; fileRef = 135F218E1E7ADBC000B9E633 /* InvariantHolder.hxx */; };
BEE45C0A2357BAB400E0AB7C /* InvariantComputation.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 135F21891E7ADBC000B9E633 /* InvariantComputation.cpp */; };
BEE45C0B2357BB9100E0AB7C /* UpdateCauchyGreenTensor.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE15AE081DC373AA00BA5D6B /* UpdateCauchyGreenTensor.cpp */; };
BEE45C0C2357BB9100E0AB7C /* UpdateCauchyGreenTensor.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE15AE091DC373AA00BA5D6B /* UpdateCauchyGreenTensor.hpp */; };
BEE45C0D2357BB9100E0AB7C /* UpdateCauchyGreenTensor.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BE15AE0A1DC373AA00BA5D6B /* UpdateCauchyGreenTensor.hxx */; };
BEE45C0E2357BB9400E0AB7C /* UpdateCauchyGreenTensor.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE15AE0E1DC373BA00BA5D6B /* UpdateCauchyGreenTensor.cpp */; };
BEE45C0F2357BB9400E0AB7C /* UpdateCauchyGreenTensor.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE15AE0F1DC373BA00BA5D6B /* UpdateCauchyGreenTensor.hpp */; };
BEE45C102357BB9400E0AB7C /* UpdateCauchyGreenTensor.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BE15AE101DC373BA00BA5D6B /* UpdateCauchyGreenTensor.hxx */; };
BEE45C112357BBC900E0AB7C /* GradientDisplacementMatrix.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 135F21801E7ADBC000B9E633 /* GradientDisplacementMatrix.cpp */; };
BEE45C122357BBD500E0AB7C /* MixedSolidIncompressibility.hxx in Headers */ = {isa = PBXBuildFile; fileRef = 37014786203B0B6700820FA4 /* MixedSolidIncompressibility.hxx */; };
BEE45C132357BBF700E0AB7C /* Helper.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE8C36A71DB50FB500919468 /* Helper.cpp */; };
BEE934791CFD8B4F00158440 /* MatrixConversion.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEE934761CFD8B4F00158440 /* MatrixConversion.cpp */; };
BEE9347A1CFD8B4F00158440 /* MatrixConversion.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BEE934771CFD8B4F00158440 /* MatrixConversion.hpp */; };
BEE9347B1CFD8B4F00158440 /* MatrixConversion.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BEE934781CFD8B4F00158440 /* MatrixConversion.hxx */; };
......@@ -11019,10 +11029,12 @@
BE7690AB22423BE3006EFE5B /* ExponentialJ1J4J6.hxx in Headers */,
BE4521471DAFBACD00807035 /* P1_to_P1b.hxx in Headers */,
BE4521441DAFBACD00807035 /* Phigher_to_P1.hxx in Headers */,
BEE45C0C2357BB9100E0AB7C /* UpdateCauchyGreenTensor.hpp in Headers */,
BE7C94711F604ED7003D2C52 /* CourtemancheRamirezNattel.hxx in Headers */,
BE45214C1DAFBACD00807035 /* P1b_to_P1.hpp in Headers */,
BEC72C422053EB6500E71849 /* QuasiIncompressibleSecondPiolaKirchhoffStressTensor.hxx in Headers */,
BEE45C052357804B00E0AB7C /* Stokes.hxx in Headers */,
BEE45C0D2357BB9100E0AB7C /* UpdateCauchyGreenTensor.hxx in Headers */,
BE8C37021DB50FB500919468 /* TransientSource.hpp in Headers */,
BE8C37191DB50FB500919468 /* AnalyticalPrestress.hxx in Headers */,
BE45213E1DAFBACD00807035 /* Phigher_to_P1.hxx in Headers */,
......@@ -11072,8 +11084,11 @@
BE7C94701F604ED7003D2C52 /* CourtemancheRamirezNattel.hpp in Headers */,
BE494103224250E800157863 /* FiberDensityJ1J4J6.hxx in Headers */,
BE8C371C1DB50FB500919468 /* InputAnalyticalPrestress.hxx in Headers */,
BEE45C102357BB9400E0AB7C /* UpdateCauchyGreenTensor.hxx in Headers */,
BEE45C122357BBD500E0AB7C /* MixedSolidIncompressibility.hxx in Headers */,
BE494102224250E800157863 /* FiberDensityJ1J4J6.hpp in Headers */,
BE8C37321DB50FB500919468 /* AnalyticalPrestress.hpp in Headers */,
BEE45C0F2357BB9400E0AB7C /* UpdateCauchyGreenTensor.hpp in Headers */,
BE8C36E81DB50FB500919468 /* Mass.hxx in Headers */,
BE4521461DAFBACD00807035 /* P1_to_P1b.hpp in Headers */,
BE7690AA22423BE3006EFE5B /* ExponentialJ1J4J6.hpp in Headers */,
......@@ -11589,6 +11604,7 @@
buildActionMask = 2147483647;
files = (
BE01DAD41E854DDD00F3EAF7 /* ExtractBlockFromGlobalVector.hpp in Headers */,
BEE45C092357BA8100E0AB7C /* InvariantHolder.hxx in Headers */,
BE01DABA1E854DB300F3EAF7 /* Helper.hpp in Headers */,
BE4521F21DAFC69D00807035 /* StVenantKirchhoff.hxx in Headers */,
BE01DAC21E854DC600F3EAF7 /* LinearLocalVariationalOperator.hpp in Headers */,
......@@ -11664,7 +11680,6 @@
BE01DAEC1E854E0100F3EAF7 /* InterpolationData.hxx in Headers */,
BE01DABF1E854DC600F3EAF7 /* BilinearLocalVariationalOperator.hxx in Headers */,
BE01DAC11E854DC600F3EAF7 /* ElementaryData.hxx in Headers */,
BE01DADB1E854DE300F3EAF7 /* InvariantHolder.hxx in Headers */,
BE01DAF81E854E1200F3EAF7 /* LocalLagrangianInterpolator.hpp in Headers */,
BE7A388D1E8D3167009DFFC3 /* GreenLagrangeTensor.hxx in Headers */,
);
......@@ -12601,15 +12616,18 @@
BE7C94721F604ED7003D2C52 /* FitzHughNagumo.cpp in Sources */,
1380F7941FCEFDF800E69537 /* GradOnGradientBasedElasticityTensor.cpp in Sources */,
BEE45C032357804B00E0AB7C /* Stokes.cpp in Sources */,
BEE45C132357BBF700E0AB7C /* Helper.cpp in Sources */,
BE8C37221DB50FB500919468 /* None.cpp in Sources */,
BEE45BFB23576FA600E0AB7C /* GradPhiGradPhi.cpp in Sources */,
BE8C37291DB50FB500919468 /* None.cpp in Sources */,
BEE45C0E2357BB9400E0AB7C /* UpdateCauchyGreenTensor.cpp in Sources */,
BE494104224250E800157863 /* FiberDensityJ1J4J6.cpp in Sources */,
13DE260B1F83D5EF00AEED6A /* MooneyRivlin.cpp in Sources */,
BE7C946F1F604ED7003D2C52 /* CourtemancheRamirezNattel.cpp in Sources */,
37014790203B15A700820FA4 /* CiarletGeymonatDeviatoric.cpp in Sources */,
37014791203B15AB00820FA4 /* LogI3Penalization.cpp in Sources */,
BE8C371E1DB50FB500919468 /* None.cpp in Sources */,
BEE45C0B2357BB9100E0AB7C /* UpdateCauchyGreenTensor.cpp in Sources */,
1380F7931FCEFDF300E69537 /* GradOnGradientBasedElasticityTensor.cpp in Sources */,
BE7690A922423BE3006EFE5B /* ExponentialJ1J4J6.cpp in Sources */,
1380F7951FCEFE9A00E69537 /* GradOnGradientBasedElasticityTensor.cpp in Sources */,
......@@ -12850,11 +12868,13 @@
BE01DAFE1E854E2500F3EAF7 /* FindCoordsOfGlobalVector.cpp in Sources */,
BE01DB041E854EB400F3EAF7 /* ExtractLocalDofValues.cpp in Sources */,
BE01DAC81E854DD600F3EAF7 /* ExtractGradientBasedBlock.cpp in Sources */,
BEE45C0A2357BAB400E0AB7C /* InvariantComputation.cpp in Sources */,
BE01DAEA1E854E0100F3EAF7 /* InterpolationData.cpp in Sources */,
BE01DAD31E854DDD00F3EAF7 /* ExtractBlockFromGlobalVector.cpp in Sources */,
BE01DAED1E854E0100F3EAF7 /* SourceOrTargetData.cpp in Sources */,
BE01DACE1E854DD600F3EAF7 /* InformationsAtQuadraturePoint.cpp in Sources */,
13B511031FCD96980080E6F6 /* DetermineExtendedUnknownList.cpp in Sources */,
BEE45C112357BBC900E0AB7C /* GradientDisplacementMatrix.cpp in Sources */,
BEF3B48B212DB18000807965 /* ForUnknownList.cpp in Sources */,
BE01DAF21E854E0D00F3EAF7 /* ComputePatternHelper.cpp in Sources */,
);
......@@ -64,10 +64,6 @@
argument = "-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/Elasticity/demo_2d.lua"
isEnabled = "NO">
</CommandLineArgument>
<CommandLineArgument
argument = "--overwrite_directory"
isEnabled = "NO">
</CommandLineArgument>
<CommandLineArgument
argument = "-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/Elasticity/demo_3d.lua"
isEnabled = "YES">
......
......@@ -68,6 +68,10 @@
argument = "-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/RivlinCube/demo_input_rivlin_cube_tetrahedra.lua"
isEnabled = "NO">
</CommandLineArgument>
<CommandLineArgument
argument = "--overwrite_directory"
isEnabled = "YES">
</CommandLineArgument>
</CommandLineArguments>
</LaunchAction>
<ProfileAction
......
......@@ -72,24 +72,24 @@ namespace MoReFEM
Internal::ConformInterpolatorNS::Local::Impl::AssertLocalNodeConsistency(p1_local_node_list, p1b_local_node_list);
#endif // NDEBUG
local_projection_matrix.Resize(static_cast<int>(p1b_local_node_list.size()) * target_data.NunknownComponent(),
static_cast<int>(p1_local_node_list.size()) * source_data.NunknownComponent());
local_projection_matrix.resize({static_cast<std::size_t>(p1b_local_node_list.size()) * target_data.NunknownComponent(),
static_cast<std::size_t>(p1_local_node_list.size()) * source_data.NunknownComponent()});
local_projection_matrix.Zero();
const auto Nnode_in_col = static_cast<int>(p1_local_node_list.size());
const auto Nnode_in_row = static_cast<int>(p1b_local_node_list.size());
const auto Nnode_in_col = p1_local_node_list.size();
const auto Nnode_in_row = p1b_local_node_list.size();
assert(Nnode_in_row >= Nnode_in_col);
// Prepare the content of the block that will be repeated for each unknown component.
LocalMatrix block(Nnode_in_row, Nnode_in_col);
block.Zero();
block.fill(0.);
for (auto i = 0; i < Nnode_in_col; ++i)
block(i, i) = 1.;
const auto Nvertex = static_cast<int>(ref_geom_elt.Nvertex());
assert(Nvertex > 0);
const double value_to_set = 1. / static_cast<double>(Nvertex);
......
......@@ -72,8 +72,8 @@ namespace MoReFEM
Internal::ConformInterpolatorNS::Local::Impl::AssertLocalNodeConsistency(p1_local_node_list, p2_local_node_list);
#endif // NDEBUG
local_projection_matrix.Resize(static_cast<int>(p2_local_node_list.size()) * target_data.NunknownComponent(),
static_cast<int>(p1_local_node_list.size()) * source_data.NunknownComponent());
local_projection_matrix.resize({p2_local_node_list.size()) * target_data.NunknownComponent(),
p1_local_node_list.size()) * source_data.NunknownComponent()});
local_projection_matrix.Zero();
......@@ -83,7 +83,7 @@ namespace MoReFEM
assert(Nnode_in_row > Nnode_in_col);
// Prepare the content of the block that will be repeated for each unknown component.
LocalMatrix block(Nnode_in_row, Nnode_in_col);
LocalMatrix block({Nnode_in_row, Nnode_in_col});
block.Zero();
for (auto i = 0; i < Nnode_in_col; ++i)
......
......@@ -68,8 +68,7 @@ namespace MoReFEM
const auto mesh_dimension = elementary_data.GetMeshDimension();
gradient_displacement_.Resize(static_cast<int>(mesh_dimension),
static_cast<int>(mesh_dimension));
gradient_displacement_.resize({ mesh_dimension, mesh_dimension });
}
......
......@@ -12,9 +12,9 @@ target_sources(${MOREFEM_OP_INSTANCES}
PRIVATE
)
#include(${CMAKE_CURRENT_LIST_DIR}/ParameterOperator/SourceList.cmake)
include(${CMAKE_CURRENT_LIST_DIR}/ParameterOperator/SourceList.cmake)
include(${CMAKE_CURRENT_LIST_DIR}/VariationalOperator/SourceList.cmake)
#include(${CMAKE_CURRENT_LIST_DIR}/HyperelasticLaws/SourceList.cmake)
include(${CMAKE_CURRENT_LIST_DIR}/HyperelasticLaws/SourceList.cmake)
#include(${CMAKE_CURRENT_LIST_DIR}/ConformInterpolator/SourceList.cmake)
#include(${CMAKE_CURRENT_LIST_DIR}/NonConformInterpolator/SourceList.cmake)
#include(${CMAKE_CURRENT_LIST_DIR}/StateOperator/SourceList.cmake)
......@@ -93,7 +93,7 @@ namespace MoReFEM
auto& elementary_data = GetNonCstElementaryData();
auto& matrix_result = elementary_data.GetNonCstMatrixResult();
matrix_result.Zero();
matrix_result.fill(0.);
const auto& unknown1_ref_felt = elementary_data.GetRefFElt(GetNthUnknown(0));
const auto& unknown2_ref_felt = elementary_data.GetRefFElt(GetNthUnknown(1));
......@@ -177,30 +177,30 @@ namespace MoReFEM
assert(dphi_test.shape(0) == Nnode_for_test_unknown1);
assert(dpsi_test.shape(0) == Nnode_for_test_unknown2);
Wrappers::Seldon::Transpose(dphi, transposed_dphi);
Wrappers::Seldon::Transpose(dpsi, transposed_dpsi);
xt::noalias(transposed_dphi) = xt::transpose(dphi);
xt::noalias(transposed_dpsi) = xt::transpose(dpsi);
block_matrix1.Zero();
block_matrix2.Zero();
block_matrix3.Zero();
block_matrix4.Zero();
block_matrix1.fill(0.);
block_matrix2.fill(0.);
block_matrix3.fill(0.);
block_matrix4.fill(0.);
Seldon::Mlt(factor1,
xt::noalias(FOO) = xt::linalg::dot((factor1,
dphi_test,
transposed_dphi,
block_matrix1);
Seldon::Mlt(factor1 + factor2,
xt::noalias(FOO) = xt::linalg::dot((factor1 + factor2,
dpsi_test,
transposed_dpsi,
block_matrix2);
Seldon::Mlt(factor1,
xt::noalias(FOO) = xt::linalg::dot((factor1,
dphi_test,
transposed_dpsi,
block_matrix3);
Seldon::Mlt(factor1,
xt::noalias(FOO) = xt::linalg::dot((factor1,
dpsi_test,
transposed_dphi,
block_matrix4);
......@@ -211,7 +211,7 @@ namespace MoReFEM
for (int component = 0 ; component < felt_space_dimension ; ++component)
norm += NumericNS::Square(tau_interpolate(component));
tau_sigma.Zero();
tau_sigma.fill(0.);
if (!(NumericNS::IsZero(norm)))
{
......@@ -219,45 +219,45 @@ namespace MoReFEM
Mlt(1. / norm, tau_sigma);
Seldon::Mlt(factor3 - factor1,
xt::noalias(FOO) = xt::linalg::dot((factor3 - factor1,
dphi_test,
tau_sigma,
dphi_test_sigma);
Seldon::MltAdd(1.,
xt::noalias(FOO) = xt::linalg::dot(Add(1.,
dphi_test_sigma,
transposed_dphi,
1.,
block_matrix1);
Seldon::Mlt(factor3 - factor1,
xt::noalias(FOO) = xt::linalg::dot((factor3 - factor1,
dphi_test,
tau_sigma,
dphi_test_sigma);
Seldon::MltAdd(1.,
xt::noalias(FOO) = xt::linalg::dot(Add(1.,
dphi_test_sigma,
transposed_dpsi,
1.,
block_matrix3);
Seldon::Mlt(factor3 - factor1,
xt::noalias(FOO) = xt::linalg::dot((factor3 - factor1,
dpsi_test,
tau_sigma,
dpsi_test_sigma);
Seldon::MltAdd(1.,
xt::noalias(FOO) = xt::linalg::dot(Add(1.,
dpsi_test_sigma,
transposed_dphi,
1.,
block_matrix4);
Seldon::Mlt(factor3 - factor1 + factor4 - factor2,
xt::noalias(FOO) = xt::linalg::dot((factor3 - factor1 + factor4 - factor2,
dpsi_test,
tau_sigma,
dpsi_test_sigma);
Seldon::MltAdd(1.,
xt::noalias(FOO) = xt::linalg::dot(Add(1.,
dpsi_test_sigma,
transposed_dpsi,
1.,
......
......@@ -165,16 +165,16 @@ namespace MoReFEM
assert(dphi.shape(0) == Nnode_for_unknown);
assert(dphi_test.shape(0) == Nnode_for_test_unknown);
Wrappers::Seldon::Transpose(dphi, transposed_dphi);
xt::noalias(transposed_dphi) = xt::transpose(dphi);
block_matrix.Zero();
Seldon::Mlt(factor1,
xt::noalias(FOO) = xt::linalg::dot((factor1,
dphi_test,
reduced_contravariant_metric_tensor,
dphi_test_contravariant_metric_tensor);
Seldon::Mlt(1.,
xt::noalias(FOO) = xt::linalg::dot((1.,
dphi_test_contravariant_metric_tensor,
transposed_dphi,
block_matrix);
......@@ -223,17 +223,17 @@ namespace MoReFEM
else
i_theta = 0.5 + std::sin(2. * theta) / (4. * theta);
Seldon::Mlt(1. - i_theta, tau_ortho_sigma);
xt::noalias(FOO) = xt::linalg::dot((1. - i_theta, tau_ortho_sigma);
Mlt(i_theta / norm, tau_sigma);
Seldon::Add(1. / norm, tau_ortho_sigma, tau_sigma);
Seldon::Mlt(factor2 - factor1,
xt::noalias(FOO) = xt::linalg::dot((factor2 - factor1,
dphi_test,
tau_sigma,
dphi_test_sigma);
Seldon::MltAdd(1.,
xt::noalias(FOO) = xt::linalg::dot(Add(1.,
dphi_test_sigma,
transposed_dphi,
1.,
......@@ -334,9 +334,9 @@ namespace MoReFEM
covariant_metric_tensor.Zero();
Wrappers::Seldon::Transpose(covariant_basis, transposed_covariant_basis);
xt::noalias(transposed_covariant_basis) = xt::transpose(covariant_basis);
Seldon::Mlt(1., transposed_covariant_basis, covariant_basis, covariant_metric_tensor);
xt::noalias(FOO) = xt::linalg::dot((1., transposed_covariant_basis, covariant_basis, covariant_metric_tensor);
double determinant = 0;
......@@ -354,7 +354,7 @@ namespace MoReFEM
contravariant_basis.Zero();
Seldon::Mlt(1., covariant_basis, contravariant_metric_tensor, contravariant_basis);
xt::noalias(FOO) = xt::linalg::dot((1., covariant_basis, contravariant_metric_tensor, contravariant_basis);
}
......
......@@ -78,7 +78,7 @@ namespace MoReFEM
auto& elementary_data = GetNonCstElementaryData();
auto& matrix_result = elementary_data.GetNonCstMatrixResult();
matrix_result.Zero();
matrix_result.fill(0.);
const auto& ref_felt = elementary_data.GetRefFElt(GetNthUnknown(0));
const auto& test_ref_felt = elementary_data.GetTestRefFElt(GetNthTestUnknown(0));
......@@ -136,11 +136,11 @@ namespace MoReFEM
assert(dphi.shape(0) == Nnode_for_unknown);
assert(dphi_test.shape(0) == Nnode_for_test_unknown);
Wrappers::Seldon::Transpose(dphi, transposed_dphi);
xt::noalias(transposed_dphi) = xt::transpose(dphi);
block_matrix.Zero();
block_matrix.fill(0.);
Seldon::Mlt(factor1,
xt::noalias(FOO) = xt::linalg::dot((factor1,
dphi_test,
transposed_dphi,
block_matrix);
......@@ -151,19 +151,19 @@ namespace MoReFEM
for (int component = 0 ; component < felt_space_dimension ; ++component)
norm += tau_interpolate(component) * tau_interpolate(component);
tau_sigma.Zero();
tau_sigma.fill(0.);
if (!(NumericNS::IsZero(norm)))
{
Wrappers::Seldon::OuterProd(tau_interpolate, tau_interpolate, tau_sigma);
Mlt(1. / norm, tau_sigma);
Seldon::Mlt(factor2 - factor1,
xt::noalias(FOO) = xt::linalg::dot((factor2 - factor1,
dphi_test,
tau_sigma,
dphi_test_sigma);
Seldon::MltAdd(1.,
xt::noalias(FOO) = xt::linalg::dot(Add(1.,
dphi_test_sigma,
transposed_dphi,
1.,
......
......@@ -120,7 +120,7 @@ namespace MoReFEM
auto& elementary_data = GetNonCstElementaryData();
auto& matrix_result = elementary_data.GetNonCstMatrixResult();
matrix_result.Zero();
matrix_result.fill(0.);
const auto& unknown1_ref_felt = elementary_data.GetRefFElt(GetNthUnknown(0));
const auto& unknown2_ref_felt = elementary_data.GetRefFElt(GetNthUnknown(1));
......@@ -221,53 +221,38 @@ namespace MoReFEM
assert(dphi_test.shape(0) == Nnode_for_test_unknown1);
assert(dpsi_test.shape(0) == Nnode_for_test_unknown2);
Wrappers::Seldon::Transpose(dphi, transposed_dphi);
Wrappers::Seldon::Transpose(dpsi, transposed_dpsi);
xt::noalias(transposed_dphi) = xt::transpose(dphi);
xt::noalias(transposed_dpsi) = xt::transpose(dpsi);
block_matrix1.Zero();
block_matrix2.Zero();
block_matrix3.Zero();
block_matrix4.Zero();
block_matrix1.fill(0.);
block_matrix2.fill(0.);
block_matrix3.fill(0.);
block_matrix4.fill(0.);
Seldon::Mlt(factor1,
dphi_test,
reduced_contravariant_metric_tensor,
dphi_test_contravariant_metric_tensor);
Seldon::Mlt(1.,
dphi_test_contravariant_metric_tensor,
transposed_dphi,
block_matrix1);
Seldon::Mlt(factor1 + factor2,
dpsi_test,
reduced_contravariant_metric_tensor,
dpsi_test_contravariant_metric_tensor);
xt::noalias(dphi_test_contravariant_metric_tensor) =
factor1 * xt::linalg::dot(dphi_test, reduced_contravariant_metric_tensor);
xt::noalias(block_matrix1) =
xt::linalg::dot(dphi_test_contravariant_metric_tensor, transposed_dphi);
Seldon::Mlt(1.,
dpsi_test_contravariant_metric_tensor,
transposed_dpsi,
block_matrix2);
xt::noalias(dpsi_test_contravariant_metric_tensor) = (factor1 + factor2) *
xt::linalg::dot(dpsi_test, reduced_contravariant_metric_tensor);
Seldon::Mlt(factor1,
dphi_test,
reduced_contravariant_metric_tensor,
dphi_test_contravariant_metric_tensor);
xt::noalias(block_matrix2) = xt::linalg::dot(dpsi_test_contravariant_metric_tensor,
transposed_dpsi);
Seldon::Mlt(1.,
dphi_test_contravariant_metric_tensor,
transposed_dpsi,
block_matrix3);
xt::noalias(dphi_test_contravariant_metric_tensor) = factor1
* xt::linalg::dot(dphi_test, reduced_contravariant_metric_tensor);
Seldon::Mlt(factor1,
dpsi_test,
reduced_contravariant_metric_tensor,
dpsi_test_contravariant_metric_tensor);
xt::noalias(block_matrix3) = xt::linalg::dot(dphi_test_contravariant_metric_tensor,
transposed_dpsi);
Seldon::Mlt(1.,
dpsi_test_contravariant_metric_tensor,
transposed_dphi,
block_matrix4);
xt::noalias(dpsi_test_contravariant_metric_tensor) =
factor1 * xt::linalg::dot(dpsi_test, reduced_contravariant_metric_tensor);
xt::noalias(block_matrix4) = xt::linalg::dot(dpsi_test_contravariant_metric_tensor, transposed_dphi);
// Fiber part of the operator.
const auto& tau_interpolate = fibers.GetValue(quad_pt, geom_elt);
......@@ -276,19 +261,19 @@ namespace MoReFEM
for (int component = 0 ; component < mesh_dimension ; ++component)
norm += tau_interpolate(component) * tau_interpolate(component);
tau_sigma.Zero();
tau_ortho_sigma.Zero();
tau_sigma.fill(0.);
tau_ortho_sigma.fill(0.);
if (!(NumericNS::IsZero(norm)))
{
// tau_interpolate_ortho
tau_interpolate_ortho.Zero();
tau_interpolate_ortho.fill(0.);
tau_interpolate_ortho(0) = tau_interpolate(1) * contravariant_basis(2,2) - tau_interpolate(2) * contravariant_basis(1,2);
tau_interpolate_ortho(1) = tau_interpolate(2) * contravariant_basis(0,2) - tau_interpolate(0) * contravariant_basis(2,2);
tau_interpolate_ortho(2) = tau_interpolate(0) * contravariant_basis(1,2) - tau_interpolate(1) * contravariant_basis(0,2);
tau_covariant_basis.Zero();
tau_ortho_covariant_basis.Zero();
tau_covariant_basis.fill(0.);
tau_ortho_covariant_basis.fill(0.);
// Projection in the covariant basis.
for (int component = 0 ; component < mesh_dimension ; ++component)
......@@ -313,50 +298,50 @@ namespace MoReFEM
else
i_theta = 0.5 + std::sin(2. * theta) / (4. * theta);
Seldon::Mlt(1. - i_theta, tau_ortho_sigma);
xt::noalias(FOO) = xt::linalg::dot((1. - i_theta, tau_ortho_sigma);
Mlt(i_theta / norm, tau_sigma);
Seldon::Add(1. / norm, tau_ortho_sigma, tau_sigma);
Seldon::Mlt((factor3 - factor1) * heterogeneous_conductivity_coefficient_at_quad,
xt::noalias(FOO) = xt::linalg::dot(((factor3 - factor1) * heterogeneous_conductivity_coefficient_at_quad,
dphi_test,
tau_sigma,
dphi_test_sigma);
Seldon::MltAdd(1.,
xt::noalias(FOO) = xt::linalg::dot(Add(1.,
dphi_test_sigma,
transposed_dphi,
1.,
block_matrix1);
Seldon::Mlt((factor3 - factor1) * heterogeneous_conductivity_coefficient_at_quad,
xt::noalias(FOO) = xt::linalg::dot(((factor3 - factor1) * heterogeneous_conductivity_coefficient_at_quad,
dphi_test,
tau_sigma,