v0.9.1
lesson4_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 
32 
33 constexpr int order = 2;
34 constexpr double young_modulus = 1;
35 constexpr double poisson_ratio = 0.25;
36 
37 #include <ElasticOps.hpp>
38 using namespace OpElasticTools;
39 
40 struct Example {
41 
42  Example(MoFEM::Interface &m_field) : mField(m_field) {}
43 
44  MoFEMErrorCode runProblem();
45 
46 private:
47  MoFEM::Interface &mField;
48 
49  MoFEMErrorCode setUP();
50  MoFEMErrorCode createCommonData();
51  MoFEMErrorCode bC();
52  MoFEMErrorCode OPs();
53  MoFEMErrorCode kspSolve();
54  MoFEMErrorCode postProcess();
55  MoFEMErrorCode checkResults();
56 
58  boost::shared_ptr<CommonData> commonDataPtr;
59 };
60 
61 //! [Run problem]
64  CHKERR setUP();
65  CHKERR createCommonData();
66  CHKERR bC();
67  CHKERR OPs();
68  CHKERR kspSolve();
69  CHKERR postProcess();
70  CHKERR checkResults();
72 }
73 //! [Run problem]
74 
75 //! [Set up problem]
78  Simple *simple = mField.getInterface<Simple>();
79  // Add field
80  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 2);
81  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 2);
82  CHKERR simple->setFieldOrder("U", order);
83  CHKERR simple->setUp();
85 }
86 //! [Set up problem]
87 
88 //! [Create common data]
91 
92  auto get_matrial_stiffens = [&](FTensor::Ddg<double, 2, 2> &t_D) {
98  t_D(i, j, k, l) = 0;
99 
100  constexpr double c = young_modulus / (1 - poisson_ratio * poisson_ratio);
101  constexpr double o = poisson_ratio * c;
102 
103  t_D(0, 0, 0, 0) = c;
104  t_D(0, 0, 1, 1) = o;
105 
106  t_D(1, 1, 0, 0) = o;
107  t_D(1, 1, 1, 1) = c;
108 
109  t_D(0, 1, 0, 1) = (1 - poisson_ratio) * c;
111  };
112 
113  commonDataPtr = boost::make_shared<CommonData>();
114  CHKERR get_matrial_stiffens(commonDataPtr->tD);
115  commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
116  commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
117  commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
119 }
120 //! [Create common data]
121 
122 //! [Boundary condition]
125 
126  auto fix_disp = [&](const std::string blockset_name) {
127  Range fix_ents;
129  if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
130  0) {
131  CHKERR mField.get_moab().get_entities_by_handle(it->meshset, fix_ents,
132  true);
133  }
134  }
135  return fix_ents;
136  };
137 
138  auto remove_ents = [&](const Range &&ents, const bool fix_x,
139  const bool fix_y) {
140  auto prb_mng = mField.getInterface<ProblemsManager>();
141  auto simple = mField.getInterface<Simple>();
143  Range verts;
144  CHKERR mField.get_moab().get_connectivity(ents, verts, true);
145  verts.merge(ents);
146  const int lo_coeff = fix_x ? 0 : 1;
147  const int hi_coeff = fix_y ? 1 : 0;
148  CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(), "U", verts,
149  lo_coeff, hi_coeff);
151  };
152 
153  CHKERR remove_ents(fix_disp("FIX_X"), true, false);
154  CHKERR remove_ents(fix_disp("FIX_Y"), false, true);
155  CHKERR remove_ents(fix_disp("FIX_ALL"), true, true);
156 
158 }
159 //! [Boundary condition]
160 
161 //! [Push operators to pipeline]
164  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
165 
166  pipeline_mng->getOpDomainLhsPipeline().push_back(
168  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetInvJacH1ForFace(invJac));
169  pipeline_mng->getOpDomainLhsPipeline().push_back(
170  new OpStiffnessMatrixLhs("U", "U", commonDataPtr));
171 
172  auto gravity = [](double x, double y) {
173  return FTensor::Tensor1<double, 2>{0., -1.};
174  };
175  pipeline_mng->getOpDomainRhsPipeline().push_back(
176  new OpForceRhs("U", commonDataPtr, gravity));
177 
178  auto integration_rule = [](int, int, int approx_order) {
179  return 2 * (approx_order - 1);
180  };
181  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
182  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
183 
185 }
186 //! [Push operators to pipeline]
187 
188 //! [Solve]
191  Simple *simple = mField.getInterface<Simple>();
192  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
193  auto solver = pipeline_mng->createKSP();
194  CHKERR KSPSetFromOptions(solver);
195  CHKERR KSPSetUp(solver);
196 
197  auto dm = simple->getDM();
198  auto D = smartCreateDMVector(dm);
199  auto F = smartVectorDuplicate(D);
200 
201  CHKERR KSPSolve(solver, F, D);
202  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
203  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
204  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
206 }
207 //! [Solve]
208 
209 //! [Postprocess results]
212  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
213 
214  pipeline_mng->getDomainLhsFE().reset();
215  auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
216  post_proc_fe->generateReferenceElementMesh();
217  post_proc_fe->getOpPtrVector().push_back(
219  post_proc_fe->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(invJac));
220  post_proc_fe->getOpPtrVector().push_back(
221  new OpCalculateVectorFieldGradient<2, 2>("U", commonDataPtr->mGradPtr));
222  post_proc_fe->getOpPtrVector().push_back(new OpStrain("U", commonDataPtr));
223  post_proc_fe->getOpPtrVector().push_back(new OpStress("U", commonDataPtr));
224  post_proc_fe->getOpPtrVector().push_back(
225  new OpPostProcElastic("U", post_proc_fe->postProcMesh,
226  post_proc_fe->mapGaussPts, commonDataPtr));
227  post_proc_fe->addFieldValuesPostProc("U");
228  pipeline_mng->getDomainRhsFE() = post_proc_fe;
229  CHKERR pipeline_mng->loopFiniteElements();
230  CHKERR post_proc_fe->writeFile("out_elastic.h5m");
232 }
233 //! [Postprocess results]
234 
235 //! [Check]
237  Simple *simple = mField.getInterface<Simple>();
238  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
240  pipeline_mng->getDomainRhsFE().reset();
241  pipeline_mng->getDomainLhsFE().reset();
242 
243  pipeline_mng->getOpDomainRhsPipeline().push_back(
245  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetInvJacH1ForFace(invJac));
246  pipeline_mng->getOpDomainRhsPipeline().push_back(
247  new OpCalculateVectorFieldGradient<2, 2>("U", commonDataPtr->mGradPtr));
248  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpStrain("U", commonDataPtr));
249  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpStress("U", commonDataPtr));
250  pipeline_mng->getOpDomainRhsPipeline().push_back(
251  new OpInternalForceRhs("U", commonDataPtr));
252 
253  auto gravity = [](double x, double y) {
254  return FTensor::Tensor1<double, 2>{0., 1.};
255  };
256  pipeline_mng->getOpDomainRhsPipeline().push_back(
257  new OpForceRhs("U", commonDataPtr, gravity));
258 
259  auto integration_rule = [](int, int, int p_data) { return 2 * (p_data - 1); };
260  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
261 
262  auto dm = simple->getDM();
263  auto res = smartCreateDMVector(dm);
264  pipeline_mng->getDomainRhsFE()->ksp_f = res;
265 
266  CHKERR VecZeroEntries(res);
267  CHKERR pipeline_mng->loopFiniteElements();
268  CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
269  CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
270  CHKERR VecAssemblyBegin(res);
271  CHKERR VecAssemblyEnd(res);
272 
273  double nrm2;
274  CHKERR VecNorm(res, NORM_2, &nrm2);
275  PetscPrintf(PETSC_COMM_WORLD, "residual = %3.4e\n", nrm2);
276  constexpr double eps = 1e-8;
277  if (nrm2 > eps)
278  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Residual is not zero");
279 
281 }
282 //! [Check]
283 
284 static char help[] = "...\n\n";
285 
286 int main(int argc, char *argv[]) {
287 
288  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
289 
290  try {
291 
292  //! [Register MoFEM discrete manager in PETSc]
293  DMType dm_name = "DMMOFEM";
294  CHKERR DMRegister_MoFEM(dm_name);
295  //! [Register MoFEM discrete manager in PETSc
296 
297  //! [Create MoAB]
298  moab::Core mb_instance; ///< mesh database
299  moab::Interface &moab = mb_instance; ///< mesh database interface
300  //! [Create MoAB]
301 
302  //! [Create MoFEM]
303  MoFEM::Core core(moab); ///< finite element database
304  MoFEM::Interface &m_field = core; ///< finite element database insterface
305  //! [Create MoFEM]
306 
307  //! [Load mesh]
308  Simple *simple = m_field.getInterface<Simple>();
309  CHKERR simple->getOptions();
310  CHKERR simple->loadFile("");
311  //! [Load mesh]
312 
313  //! [Example]
314  Example ex(m_field);
315  CHKERR ex.runProblem();
316  //! [Example]
317  }
318  CATCH_ERRORS;
319 
321 }
boost::shared_ptr< FEMethod > & getDomainLhsFE()
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Deprecated interface functions.
Problem manager is used to build and partition problems.
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
int main(int argc, char *argv[])
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:445
static constexpr int approx_order
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
MatrixDouble invJac
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
Example(MoFEM::Interface &m_field)
Core (interface) class.
Definition: Core.hpp:50
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
Definition: AuxPETSc.hpp:238
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
constexpr double young_modulus
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
constexpr double poisson_ratio
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:939
MoFEMErrorCode runProblem()
[Operator]
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEMErrorCode OPs()
[Boundary condition]
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
static char help[]
[Check]
MoFEMErrorCode checkResults()
[Print results]
Calculate inverse of jacobian for face element.
MoFEMErrorCode bC()
[Create common data]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
boost::shared_ptr< FEMethod > & getDomainRhsFE()
DataForcesAndSourcesCore::EntData EntData
PipelineManager interface.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
#define CHKERR
Inline error check.
Definition: definitions.h:602
MoFEMErrorCode setUP()
[Set up problem]
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Transform local reference derivatives of shape functions to global derivatives.
constexpr int order
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
MoFEMErrorCode postProcess()
[Integrate]
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: ElasticOps.hpp:86
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
continuous field
Definition: definitions.h:177
static const double eps
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
FTensor::Index< 'l', 2 > l
Definition: ContactOps.hpp:29
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:67
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26