Commit 684f7873 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1292 Stokes now works with Xtensor (with the Stokes operator only at the moment).

parent ec417ee9
......@@ -1468,6 +1468,12 @@
BEE45C0023576FD000E0AB7C /* ThreeDimensionalInitialCondition.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEACCC2B1E11D05600CBA4F2 /* ThreeDimensionalInitialCondition.cpp */; };
BEE45C0123576FD000E0AB7C /* ThreeDimensionalInitialCondition.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BEACCC2C1E11D05600CBA4F2 /* ThreeDimensionalInitialCondition.hpp */; };
BEE45C0223576FD000E0AB7C /* ThreeDimensionalInitialCondition.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BEACCC2D1E11D05600CBA4F2 /* ThreeDimensionalInitialCondition.hxx */; };
BEE45C032357804B00E0AB7C /* Stokes.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE8C36721DB50FB500919468 /* Stokes.cpp */; };
BEE45C042357804B00E0AB7C /* Stokes.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE8C36731DB50FB500919468 /* Stokes.hpp */; };
BEE45C052357804B00E0AB7C /* Stokes.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BE8C36741DB50FB500919468 /* Stokes.hxx */; };
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 */; };
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 */; };
......@@ -11016,6 +11022,7 @@
BE7C94711F604ED7003D2C52 /* CourtemancheRamirezNattel.hxx in Headers */,
BE45214C1DAFBACD00807035 /* P1b_to_P1.hpp in Headers */,
BEC72C422053EB6500E71849 /* QuasiIncompressibleSecondPiolaKirchhoffStressTensor.hxx in Headers */,
BEE45C052357804B00E0AB7C /* Stokes.hxx in Headers */,
BE8C37021DB50FB500919468 /* TransientSource.hpp in Headers */,
BE8C37191DB50FB500919468 /* AnalyticalPrestress.hxx in Headers */,
BE45213E1DAFBACD00807035 /* Phigher_to_P1.hxx in Headers */,
......@@ -11033,10 +11040,12 @@
BE4521491DAFBACD00807035 /* P1_to_P2.hpp in Headers */,
BE8C36E71DB50FB500919468 /* Mass.hpp in Headers */,
BE8C37371DB50FB500919468 /* Hyperelasticity.hpp in Headers */,
BEE45C042357804B00E0AB7C /* Stokes.hpp in Headers */,
BE8C37181DB50FB500919468 /* AnalyticalPrestress.hpp in Headers */,
BE8C37271DB50FB500919468 /* PartialSpecialization.hpp in Headers */,
BE8C36E01DB50FB500919468 /* GradOnGradientBasedElasticityTensor.hxx in Headers */,
BE8C370A1DB50FB500919468 /* FollowingPressure.hpp in Headers */,
BEE45C072357804F00E0AB7C /* Stokes.hpp in Headers */,
BE8C372A1DB50FB500919468 /* None.hpp in Headers */,
BE8C37391DB50FB500919468 /* None.hpp in Headers */,
BE8C37251DB50FB500919468 /* Helper.hpp in Headers */,
......@@ -11049,6 +11058,7 @@
BE8C37381DB50FB500919468 /* None.hpp in Headers */,
BE8C37211DB50FB500919468 /* Hyperelasticity.hxx in Headers */,
BE8C371F1DB50FB500919468 /* None.hpp in Headers */,
BEE45C082357804F00E0AB7C /* Stokes.hxx in Headers */,
BE8C37041DB50FB500919468 /* TransientSource.hpp in Headers */,
BE8C36F51DB50FB500919468 /* Mass.hpp in Headers */,
BE8C371B1DB50FB500919468 /* InputAnalyticalPrestress.hpp in Headers */,
......@@ -12590,6 +12600,7 @@
BE45213F1DAFBACD00807035 /* Check.cpp in Sources */,
BE7C94721F604ED7003D2C52 /* FitzHughNagumo.cpp in Sources */,
1380F7941FCEFDF800E69537 /* GradOnGradientBasedElasticityTensor.cpp in Sources */,
BEE45C032357804B00E0AB7C /* Stokes.cpp in Sources */,
BE8C37221DB50FB500919468 /* None.cpp in Sources */,
BEE45BFB23576FA600E0AB7C /* GradPhiGradPhi.cpp in Sources */,
BE8C37291DB50FB500919468 /* None.cpp in Sources */,
......@@ -12605,6 +12616,7 @@
BE8C36E61DB50FB500919468 /* Mass.cpp in Sources */,
135F59FD1F851EB300655154 /* ExponentialJ1J4.cpp in Sources */,
BE8C36F41DB50FB500919468 /* Mass.cpp in Sources */,
BEE45C062357804F00E0AB7C /* Stokes.cpp in Sources */,
BEE45BFA23576FA100E0AB7C /* GradPhiGradPhi.cpp in Sources */,
);
runOnlyForDeploymentPostprocessing = 0;
......@@ -64,6 +64,10 @@
argument = "-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/Stokes/demo.lua"
isEnabled = "YES">
</CommandLineArgument>
<CommandLineArgument
argument = "--overwrite_directory"
isEnabled = "YES">
</CommandLineArgument>
</CommandLineArguments>
</LaunchAction>
<ProfileAction
......
......@@ -94,7 +94,7 @@ namespace MoReFEM
auto& elementary_data = GetNonCstElementaryData();
auto& matrix_result = elementary_data.GetNonCstMatrixResult();
matrix_result.Zero();
matrix_result.fill(0.);
const auto& scalar_ref_felt =
elementary_data.GetRefFElt(GetNthUnknown(EnumUnderlyingType(UnknownIndex::scalar)));
......@@ -106,13 +106,13 @@ namespace MoReFEM
const auto& vectorial_test_ref_felt =
elementary_data.GetTestRefFElt(GetNthTestUnknown(EnumUnderlyingType(TestUnknownIndex::vectorial)));
const auto Nnode_velocity = static_cast<int>(vectorial_ref_felt.Nnode());
const auto Ndof_velocity = static_cast<int>(vectorial_ref_felt.Ndof());
const auto Nnode_pressure = static_cast<int>(scalar_ref_felt.Nnode());
const auto Nnode_velocity = static_cast<std::size_t>(vectorial_ref_felt.Nnode());
const auto Ndof_velocity = static_cast<std::size_t>(vectorial_ref_felt.Ndof());
const auto Nnode_pressure = static_cast<std::size_t>(scalar_ref_felt.Nnode());
const auto Nnode_test_velocity = static_cast<int>(vectorial_test_ref_felt.Nnode());
const auto Nnode_test_pressure = static_cast<int>(scalar_test_ref_felt.Nnode());
const auto Ndof_test_velocity = static_cast<int>(vectorial_test_ref_felt.Ndof());
const auto Nnode_test_velocity = static_cast<std::size_t>(vectorial_test_ref_felt.Nnode());
const auto Nnode_test_pressure = static_cast<std::size_t>(scalar_test_ref_felt.Nnode());
const auto Ndof_test_velocity = static_cast<std::size_t>(vectorial_test_ref_felt.Ndof());
const auto& infos_at_quad_pt_list =
elementary_data.GetInformationsAtQuadraturePointList();
......@@ -123,9 +123,9 @@ namespace MoReFEM
const auto& fluid_viscosity = GetFluidViscosity();
const auto& geom_elt = elementary_data.GetCurrentGeomElt();
const auto Ncomponent = static_cast<int>(vectorial_ref_felt.Ncomponent());
const auto Ncomponent = static_cast<std::size_t>(vectorial_ref_felt.Ncomponent());
assert(Ncomponent == static_cast<int>(vectorial_test_ref_felt.Ncomponent()));
assert(Ncomponent == static_cast<std::size_t>(vectorial_test_ref_felt.Ncomponent()));
for (const auto& infos_at_quad_pt : infos_at_quad_pt_list)
......@@ -154,22 +154,19 @@ namespace MoReFEM
scalar_test_ref_felt); // only on p*
// Part in (v* p)
for (int node_test_velocity_index = 0; node_test_velocity_index < Nnode_test_velocity; ++node_test_velocity_index)
for (auto node_test_velocity_index = 0ul; node_test_velocity_index < Nnode_test_velocity; ++node_test_velocity_index)
{
for (int node_pressure_index = 0; node_pressure_index < Nnode_pressure; ++node_pressure_index)
for (auto node_pressure_index = 0ul; node_pressure_index < Nnode_pressure; ++node_pressure_index)
{
// FIX THAT! int dof_test_velocity_index = node_test_velocity_index;
// A test must be written for scalar div vectorial...
const int& dof_pressure = node_pressure_index; // alias
const auto dof_pressure = node_pressure_index; // alias
// The two terms below give the same result but in a different manner. It is just here to show how to use them.
//const auto test_pressure_term = test_felt_phi(Nnode_test_velocity + node_test_pressure_index);
const auto pressure_term = pressure_felt_phi(node_pressure_index);
for (int component = 0; component < Ncomponent; ++component)
for (auto component = 0ul; component < Ncomponent; ++component)
{
int dof_test_velocity_index = node_test_velocity_index + Nnode_test_velocity * component;
auto dof_test_velocity_index = node_test_velocity_index + Nnode_test_velocity * component;
const double product =
factor * pressure_term * test_gradient_felt_phi(node_test_velocity_index, component);
......@@ -182,19 +179,19 @@ namespace MoReFEM
}
// Part in (p* v)
for (int node_test_pressure_index = 0; node_test_pressure_index < Nnode_test_pressure; ++node_test_pressure_index)
for (auto node_test_pressure_index = 0ul; node_test_pressure_index < Nnode_test_pressure; ++node_test_pressure_index)
{
const int& dof_test_pressure = node_test_pressure_index; // alias
const auto dof_test_pressure = node_test_pressure_index;
// The two terms below give the same result but in a different manner. It is just here to show how to use them.
// const auto pressure_term = felt_phi(Nnode_velocity + node_pressure_index);
const auto test_pressure_term = test_pressure_felt_phi(node_test_pressure_index);
for (int node_velocity_index = 0; node_velocity_index < Nnode_velocity; ++node_velocity_index)
for (auto node_velocity_index = 0ul; node_velocity_index < Nnode_velocity; ++node_velocity_index)
{
int dof_velocity_index = node_velocity_index;
auto dof_velocity_index = node_velocity_index;
for (int component = 0; component < Ncomponent; ++component, dof_velocity_index += Nnode_velocity)
for (auto component = 0ul; component < Ncomponent; ++component, dof_velocity_index += Nnode_velocity)
{
const double product =
factor * test_pressure_term * gradient_felt_phi(node_velocity_index, component);
......@@ -216,7 +213,7 @@ namespace MoReFEM
const auto& dphi_test_velocity = ExtractSubMatrix(test_unknown_data.GetGradientFEltPhi(),
vectorial_test_ref_felt);
velocity_block_matrix.Zero();
velocity_block_matrix.fill(0.);
// First compute the content of the block matrix.
assert(dphi_velocity.shape(0) == Nnode_velocity);
......@@ -225,21 +222,19 @@ namespace MoReFEM
assert(dphi_test_velocity.shape(0) == Nnode_test_velocity);
assert(dphi_test_velocity.shape(1) == Ncomponent);
Wrappers::Seldon::Transpose(dphi_velocity, transposed_dphi_velocity);
xt::noalias(transposed_dphi_velocity) = xt::transpose(dphi_velocity);
Seldon::Mlt(factor,
dphi_test_velocity,
transposed_dphi_velocity,
velocity_block_matrix);
xt::noalias(velocity_block_matrix) = factor * xt::linalg::dot(dphi_test_velocity,
transposed_dphi_velocity);
// Then report it into the elementary matrix.
for (int m = 0; m < Nnode_test_velocity; ++m)
for (auto m = 0ul; m < Nnode_test_velocity; ++m)
{
for (int n = 0; n < Nnode_velocity; ++n)
for (auto n = 0ul; n < Nnode_velocity; ++n)
{
const double value = velocity_block_matrix(m, n);
for (int component = 0; component < Ncomponent; ++component)
for (auto component = 0ul; component < Ncomponent; ++component)
{
const auto shift_row = Nnode_test_velocity * component;
const auto shift_col = Nnode_velocity * component;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment