v0.13.0
elastic.cpp
Go to the documentation of this file.
1 /**
2  * \file lesson4_elastic.cpp
3  * \example lesson4_elastic.cpp
4  *
5  * Plane stress elastic problem
6  *
7  */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #include <MoFEM.hpp>
24 
25 using namespace MoFEM;
26 
27 template <int DIM> struct ElementsAndOps {};
28 
29 template <> struct ElementsAndOps<2> {
35 };
36 
37 template <> struct ElementsAndOps<3> {
43 };
44 
45 //! [Define dimension]
46 constexpr int SPACE_DIM = 2; //< Space dimension of problem, mesh
47 //! [Define dimension]
48 
55 
57  GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
59  GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 0>;
62 
63 constexpr double young_modulus = 100;
64 constexpr double poisson_ratio = 0.3;
65 constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
66 constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
67 
68 #include <OpPostProcElastic.hpp>
69 using namespace Tutorial;
70 
71 struct Example {
72 
73  Example(MoFEM::Interface &m_field) : mField(m_field) {}
74 
76 
77 private:
78  MoFEM::Interface &mField;
79 
88 
89  boost::shared_ptr<MatrixDouble> matGradPtr;
90  boost::shared_ptr<MatrixDouble> matStrainPtr;
91  boost::shared_ptr<MatrixDouble> matStressPtr;
92  boost::shared_ptr<MatrixDouble> matDPtr;
93  boost::shared_ptr<MatrixDouble> bodyForceMatPtr;
94 };
95 
96 //! [Create common data]
99 
100  //! [Calculate elasticity tensor]
101  auto set_material_stiffness = [&]() {
108  constexpr double A =
109  (SPACE_DIM == 2) ? 2 * shear_modulus_G /
110  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
111  : 1;
112  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
113  t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
114  A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
115  t_kd(i, j) * t_kd(k, l);
117  };
118  //! [Calculate elasticity tensor]
119 
120  //! [Define gravity vector]
121  auto set_body_force = [&]() {
124  auto t_force = getFTensor1FromMat<SPACE_DIM, 0>(*bodyForceMatPtr);
125  t_force(i) = 0;
126  t_force(1) = -1;
128  };
129  //! [Define gravity vector]
130 
131  //! [Initialise containers for commonData]
132  matGradPtr = boost::make_shared<MatrixDouble>();
133  matStrainPtr = boost::make_shared<MatrixDouble>();
134  matStressPtr = boost::make_shared<MatrixDouble>();
135  matDPtr = boost::make_shared<MatrixDouble>();
136  bodyForceMatPtr = boost::make_shared<MatrixDouble>();
137 
138  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
139  matDPtr->resize(size_symm * size_symm, 1);
140  bodyForceMatPtr->resize(SPACE_DIM, 1);
141  //! [Initialise containers for commonData]
142 
143  CHKERR set_material_stiffness();
144  CHKERR set_body_force();
145 
147 }
148 //! [Create common data]
149 
150 //! [Run problem]
153  CHKERR readMesh();
154  CHKERR setupProblem();
155  CHKERR createCommonData();
156  CHKERR boundaryCondition();
157  CHKERR assembleSystem();
158  CHKERR solveSystem();
159  CHKERR outputResults();
160  CHKERR checkResults();
162 }
163 //! [Run problem]
164 
165 //! [Read mesh]
168  auto simple = mField.getInterface<Simple>();
169  CHKERR simple->getOptions();
170  CHKERR simple->loadFile();
172 }
173 //! [Read mesh]
174 
175 //! [Set up problem]
178  Simple *simple = mField.getInterface<Simple>();
179  // Add field
180  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
181  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
182  int order = 3;
183  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
184  CHKERR simple->setFieldOrder("U", order);
185  CHKERR simple->setUp();
187 }
188 //! [Set up problem]
189 
190 //! [Boundary condition]
193  auto *pipeline_mng = mField.getInterface<PipelineManager>();
194  auto simple = mField.getInterface<Simple>();
195  auto bc_mng = mField.getInterface<BcManager>();
196 
197  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
198  "U", 0, 0);
199  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
200  "U", 1, 1);
201  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
202  "U", 2, 2);
203  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
204  "U", 0, 3);
205 
206  //! [Pushing gravity load operator]
207  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpBodyForce(
208  "U", bodyForceMatPtr, [](double, double, double) { return 1.; }));
209  //! [Pushing gravity load operator]
210 
212 }
213 //! [Boundary condition]
214 
215 //! [Push operators to pipeline]
218  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
219 
220  if (SPACE_DIM == 2) {
221  auto det_ptr = boost::make_shared<VectorDouble>();
222  auto jac_ptr = boost::make_shared<MatrixDouble>();
223  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
224  pipeline_mng->getOpDomainLhsPipeline().push_back(
225  new OpCalculateHOJacForFace(jac_ptr));
226  pipeline_mng->getOpDomainLhsPipeline().push_back(
227  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
228  pipeline_mng->getOpDomainLhsPipeline().push_back(
229  new OpSetInvJacH1ForFace(inv_jac_ptr));
230  }
231  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpK("U", "U", matDPtr));
232 
233  auto integration_rule = [](int, int, int approx_order) {
234  return 2 * (approx_order - 1);
235  };
238 
240 }
241 //! [Push operators to pipeline]
242 
243 //! [Solve]
246  auto *simple = mField.getInterface<Simple>();
247  auto *pipeline_mng = mField.getInterface<PipelineManager>();
248  auto solver = pipeline_mng->createKSP();
249  CHKERR KSPSetFromOptions(solver);
250  CHKERR KSPSetUp(solver);
251 
252  auto dm = simple->getDM();
253  auto D = smartCreateDMVector(dm);
254  auto F = smartVectorDuplicate(D);
255 
256  CHKERR KSPSolve(solver, F, D);
257  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
258  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
259  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
261 }
262 //! [Solve]
263 
264 //! [Postprocess results]
267  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
268  auto det_ptr = boost::make_shared<VectorDouble>();
269  auto jac_ptr = boost::make_shared<MatrixDouble>();
270  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
271  pipeline_mng->getDomainLhsFE().reset();
272  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
273  post_proc_fe->generateReferenceElementMesh();
274  if (SPACE_DIM == 2) {
275  post_proc_fe->getOpPtrVector().push_back(
276  new OpCalculateHOJacForFace(jac_ptr));
277  post_proc_fe->getOpPtrVector().push_back(
278  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
279  post_proc_fe->getOpPtrVector().push_back(
280  new OpSetInvJacH1ForFace(inv_jac_ptr));
281  }
282  post_proc_fe->getOpPtrVector().push_back(
284  matGradPtr));
285  post_proc_fe->getOpPtrVector().push_back(
286  new OpSymmetrizeTensor<SPACE_DIM>("U", matGradPtr, matStrainPtr));
287  post_proc_fe->getOpPtrVector().push_back(
289  "U", matStrainPtr, matStressPtr, matDPtr));
290  post_proc_fe->getOpPtrVector().push_back(new OpPostProcElastic<SPACE_DIM>(
291  "U", post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts, matStrainPtr,
292  matStressPtr));
293  post_proc_fe->addFieldValuesPostProc("U");
294  pipeline_mng->getDomainRhsFE() = post_proc_fe;
295  CHKERR pipeline_mng->loopFiniteElements();
296  CHKERR post_proc_fe->writeFile("out_elastic.h5m");
298 }
299 //! [Postprocess results]
300 
301 //! [Check]
303  MOFEM_LOG_CHANNEL("WORLD");
304  Simple *simple = mField.getInterface<Simple>();
305  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
307  pipeline_mng->getDomainRhsFE().reset();
308  pipeline_mng->getDomainLhsFE().reset();
309 
310  if (SPACE_DIM == 2) {
311  auto det_ptr = boost::make_shared<VectorDouble>();
312  auto jac_ptr = boost::make_shared<MatrixDouble>();
313  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
314  pipeline_mng->getOpDomainRhsPipeline().push_back(
315  new OpCalculateHOJacForFace(jac_ptr));
316  pipeline_mng->getOpDomainRhsPipeline().push_back(
317  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
318  pipeline_mng->getOpDomainRhsPipeline().push_back(
319  new OpSetInvJacH1ForFace(inv_jac_ptr));
320  }
321  pipeline_mng->getOpDomainRhsPipeline().push_back(
323  matGradPtr));
324  pipeline_mng->getOpDomainRhsPipeline().push_back(
325  new OpSymmetrizeTensor<SPACE_DIM>("U", matGradPtr, matStrainPtr));
326  pipeline_mng->getOpDomainRhsPipeline().push_back(
328  "U", matStrainPtr, matStressPtr, matDPtr));
329  pipeline_mng->getOpDomainRhsPipeline().push_back(
330  new OpInternalForce("U", matStressPtr));
331  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpBodyForce(
332  "U", bodyForceMatPtr, [](double, double, double) { return -1.; }));
333 
334  auto integration_rule = [](int, int, int p_data) { return 2 * (p_data - 1); };
335  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
336 
337  auto dm = simple->getDM();
338  auto res = smartCreateDMVector(dm);
339  pipeline_mng->getDomainRhsFE()->ksp_f = res;
340 
341  CHKERR VecZeroEntries(res);
342  CHKERR pipeline_mng->loopFiniteElements();
343  CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
344  CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
345  CHKERR VecAssemblyBegin(res);
346  CHKERR VecAssemblyEnd(res);
347 
348  double nrm2;
349  CHKERR VecNorm(res, NORM_2, &nrm2);
350  MOFEM_LOG_C("WORLD", Sev::verbose, "residual = %3.4e\n", nrm2);
351  constexpr double eps = 1e-8;
352  if (nrm2 > eps)
353  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Residual is not zero");
354 
356 }
357 //! [Check]
358 
359 static char help[] = "...\n\n";
360 
361 int main(int argc, char *argv[]) {
362 
363  // Initialisation of MoFEM/PETSc and MOAB data structures
364  const char param_file[] = "param_file.petsc";
365  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
366 
367  try {
368 
369  //! [Register MoFEM discrete manager in PETSc]
370  DMType dm_name = "DMMOFEM";
371  CHKERR DMRegister_MoFEM(dm_name);
372  //! [Register MoFEM discrete manager in PETSc
373 
374  //! [Create MoAB]
375  moab::Core mb_instance; ///< mesh database
376  moab::Interface &moab = mb_instance; ///< mesh database interface
377  //! [Create MoAB]
378 
379  //! [Create MoFEM]
380  MoFEM::Core core(moab); ///< finite element database
381  MoFEM::Interface &m_field = core; ///< finite element database insterface
382  //! [Create MoFEM]
383 
384  //! [Example]
385  Example ex(m_field);
386  CHKERR ex.runProblem();
387  //! [Example]
388  }
389  CATCH_ERRORS;
390 
392 }
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static const double eps
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
int main(int argc, char *argv[])
Definition: elastic.cpp:361
EntitiesFieldData::EntData EntData
[Define dimension]
Definition: elastic.cpp:49
static char help[]
[Check]
Definition: elastic.cpp:359
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBodyForce
Definition: elastic.cpp:59
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:46
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
Definition: elastic.cpp:61
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
Definition: elastic.cpp:57
constexpr double poisson_ratio
Definition: elastic.cpp:64
constexpr double shear_modulus_G
Definition: elastic.cpp:66
constexpr double bulk_modulus_K
Definition: elastic.cpp:65
constexpr double young_modulus
Definition: elastic.cpp:63
auto integration_rule
constexpr auto t_kd
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:287
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
double A
static constexpr int approx_order
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
PipelineManager::FaceEle DomainEle
[Example]
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
boost::shared_ptr< MatrixDouble > bodyForceMatPtr
Definition: elastic.cpp:93
Example(MoFEM::Interface &m_field)
Definition: elastic.cpp:73
MoFEMErrorCode runProblem()
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Postprocess on face.
Post processing.
[Class definition]