v0.14.0
mwls_at_gauss_point.cpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 #ifdef WITH_TETGEN
16 
17 #include <tetgen.h>
18 #ifdef REAL
19 #undef REAL
20 #endif
21 
22 #endif
23 
24 #include <moab/SpatialLocator.hpp>
25 #include <MoFEM.hpp>
26 using namespace MoFEM;
27 #include <BasicFiniteElements.hpp>
28 #include <Mortar.hpp>
29 #include <NeoHookean.hpp>
30 #include <Hooke.hpp>
31 
32 #define BOOST_UBLAS_NDEBUG 1
33 
34 #include <CrackFrontElement.hpp>
35 #include <ComplexConstArea.hpp>
36 #include <MWLS.hpp>
37 #include <GriffithForceElement.hpp>
38 #include <VolumeLengthQuality.hpp>
39 #include <CrackPropagation.hpp>
40 
41 
42 static char help[] = "Testing moving weighted least square approximation"
43  "n\n";
44 
45 #define MWLSLog \
46  MOFEM_LOG_CHANNEL("WORLD"); \
47  MOFEM_LOG_FUNCTION(); \
48  MOFEM_LOG_TAG("WORLD", "mwls_at_gauss_point");
49 
50 using namespace FractureMechanics;
51 
52 bool CrackPropagation::debug = false;
53 
54 bool CrackPropagation::parallelReadAndBroadcast =
55  false; ///< That is unstable development, for
56  ///< some meshses (propably generated by
57  ///< cubit) this is not working. Error can
58  ///< be attributed to bug in MOAB.
59 
60 int main(int argc, char *argv[]) {
61  const char param_file[] = "param_file.petsc";
62  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
63 
64  try {
65 
66  int order = 1;
67  PetscBool test;
68  double young_modulus = 18000;
69  double poisson_ratio = 0.2;
70 
71  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "mwls_at_gauss_point",
72  "none");
73  CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
74  PETSC_NULL);
75  CHKERR PetscOptionsBool("-test", "test code", "", test, &test, PETSC_NULL);
76  CHKERR PetscOptionsScalar("-young_modulus",
77  "fix nodes below given threshold", "",
78  young_modulus, &young_modulus, PETSC_NULL);
79  CHKERR PetscOptionsScalar("-poisson_ratio",
80  "fix nodes below given threshold", "",
81  poisson_ratio, &poisson_ratio, PETSC_NULL);
82 
83  ierr = PetscOptionsEnd();
84  CHKERRG(ierr);
85 
86  // Register DM Manager
87  DMType dm_name = "DMMOFEM";
88  CHKERR DMRegister_MoFEM(dm_name);
89 
90  moab::Core mb_instance;
91  moab::Interface &moab = mb_instance;
92 
93  MoFEM::Core core(moab);
94  MoFEM::Interface &m_field = core;
95 
96  // *** Create simple interface to get data from element Gauss point
97  Simple *simple_interface_ptr;
98  CHKERR m_field.getInterface(simple_interface_ptr);
99  // Get options from command line
100  CHKERR simple_interface_ptr->getOptions();
101 
102  // Load file
103  CHKERR simple_interface_ptr->loadFile("", "");
104  Range ents;
105  EntityHandle whole_mesh = 0;
106  BitRefLevel bitLevel;
107  bitLevel.set(0);
108  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(0, 3,
109  bitLevel);
110 
111  CHKERR m_field.get_moab().get_entities_by_dimension(whole_mesh, 3, ents,
112  true);
113  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(ents, bitLevel,
114  false);
115 
116  // *** Add fields
117  CHKERR simple_interface_ptr->addDomainField("MESH_NODE_POSITIONS", H1,
118  AINSWORTH_LEGENDRE_BASE, 3);
119  CHKERR simple_interface_ptr->addDomainField("SPATIAL_POSITION", H1,
120  AINSWORTH_LEGENDRE_BASE, 3);
121  // *** Set field order
122  CHKERR simple_interface_ptr->setFieldOrder("MESH_NODE_POSITIONS", order);
123  CHKERR simple_interface_ptr->setFieldOrder("SPATIAL_POSITION", order);
124  CHKERR simple_interface_ptr->buildFields();
125 
126  Projection10NodeCoordsOnField ent_method_spatial1(m_field,
127  "SPATIAL_POSITION");
128  CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method_spatial1, NOISY);
129 
130  Projection10NodeCoordsOnField ent_method_material2(m_field,
131  "MESH_NODE_POSITIONS");
132 
133  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material2, 0);
134  CHKERR simple_interface_ptr->setUp(PETSC_FALSE);
135 
136  // *** Create elastic element data structure
137  auto elastic_fe =
138  boost::make_shared<NonlinearElasticElement>(m_field, ELASTIC_TAG);
139 
140  {
141 
142  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
143  m_field, BLOCKSET | MAT_ELASTICSET, it)) {
144  // Get data from block
145  Mat_Elastic mydata;
146  CHKERR it->getAttributeDataStructure(mydata);
147  int id = it->getMeshsetId();
148  EntityHandle meshset = it->getMeshset();
149  CHKERR m_field.get_moab().get_entities_by_type(
150  meshset, MBTET, elastic_fe->setOfBlocks[id].tEts, true);
151  elastic_fe->setOfBlocks[id].iD = id;
152 
153  elastic_fe->setOfBlocks[id].materialDoublePtr =
154  boost::make_shared<Hooke<double>>();
155  elastic_fe->setOfBlocks[id].materialAdoublePtr =
156  boost::make_shared<Hooke<adouble>>();
157  }
158 
159  elastic_fe->commonData.spatialPositions = "SPATIAL_POSITION";
160  elastic_fe->commonData.meshPositions = "MESH_NODE_POSITIONS";
161 
162  // *** Create instance for moving least square approximation
163 
164  CrackPropagation cp(m_field, 3, 2);
165  CHKERR cp.getOptions();
166  auto mwls_approx = boost::make_shared<MWLSApprox>(m_field);
167 
168  {
169 
170  CHKERR mwls_approx->loadMWLSMesh(cp.mwlsApproxFile.c_str());
171 
172  MeshsetsManager *block_manager_ptr;
173  CHKERR m_field.getInterface(block_manager_ptr);
174  EntityHandle approx_meshset = 0;
175  if (cp.residualStressBlock != -1) {
176  CHKERR block_manager_ptr->getMeshset(cp.residualStressBlock, BLOCKSET,
177  approx_meshset);
178  }
179  // get tets on which stress field is going to be approximated
180  Range tets;
181  CHKERR m_field.get_moab().get_entities_by_type(whole_mesh, MBTET, tets,
182  true);
183  EntityHandle meshset;
184  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, meshset);
185  CHKERR m_field.get_moab().add_entities(meshset, tets);
186 
187  // Create tags to approximate data
188  std::vector<Tag> tag_handles;
189  // Get handles for all tags defined in the mesh instance.
190  CHKERR mwls_approx->mwlsMoab.tag_get_tags(tag_handles);
191  ublas::vector<Tag> tag_approx_handles(tag_handles.size());
192  ublas::matrix<Tag> tag_approx_handles_diff(3, tag_handles.size());
193  ublas::matrix<Tag> tag_approx_handles_diff_diff(9, tag_handles.size());
194  double def_vals[9];
195  fill(def_vals, &def_vals[9], 0);
196 
197  ublas::vector<Tag> tag_difference_handles(tag_handles.size());
198  double def_differece_vals[9];
199  fill(def_differece_vals, &def_differece_vals[9], 0);
200 
201  ublas::vector<Tag> tag_error_handles(tag_handles.size());
202  double def_error_vals = 0;
203 
204  ublas::vector<Tag> tag_hydro_static_error_handles(tag_handles.size());
205  double def_hydro_static_error_vals = 0;
206 
207  ublas::vector<Tag> tag_deviatoric_difference_handles(
208  tag_handles.size());
209  double def_deviatoric_differece_vals[9];
210  fill(def_deviatoric_differece_vals, &def_deviatoric_differece_vals[9],
211  0);
212 
213  ublas::vector<Tag> tag_deviatoric_error_handles(tag_handles.size());
214  double def_deviatoric_error_vals = 0;
215 
216  // interate over all tags
217  for (int t = 0; t != tag_handles.size(); t++) {
218  // get tag name
219  std::string name;
220  CHKERR mwls_approx->mwlsMoab.tag_get_name(tag_handles[t], name);
221  // check if tag is the given tag
222  if (name.compare(cp.mwlsStressTagName) == 0 &&
223  name.compare("RHO") != 0) {
224 
225  // create elements instances
226  Range ref_nodes; // = crackFrontNodes;
227  Range ref_node_edges; // = crackFrontNodesEdges;
228  Range ref_element; // = crackFrontElements
229  PetscBool boolean = PETSC_FALSE;
230 
231  auto error_volume_fe = boost::shared_ptr<CrackFrontElement>(
232  new CrackFrontElement(m_field, boolean, ref_nodes,
233  ref_node_edges, ref_element, boolean));
234 
235  {
236 
237  error_volume_fe->meshPositionsFieldName = "NONE";
238  error_volume_fe->addToRule = 0;
239 
240  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(
241  new MatrixDouble());
242  error_volume_fe->getOpPtrVector().push_back(
243  new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS",
244  mat_pos_at_pts_ptr));
245 
246  boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr(
247  new MatrixDouble());
248 
249  error_volume_fe->getOpPtrVector().push_back(
250  new MWLSApprox::OpMWLSStressAtGaussPts(
251  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
252  error_volume_fe, mwls_approx, cp.mwlsStressTagName, false,
253  false));
254 
255  boost::shared_ptr<VectorDouble> vec_gp_weights(
256  new VectorDouble());
257 
258  error_volume_fe->getOpPtrVector().push_back(
259  new MWLSApprox::OpMWLSStressAndErrorsAtGaussPts(
260  mwls_approx, cp.mwlsStressTagName, m_field));
261 
262  CHKERR block_manager_ptr->getEntitiesByDimension(
263  cp.residualStressBlock, BLOCKSET, 3, mwls_approx->tetsInBlock,
264  true);
265  // No member name getEntitiesByDimension in Simple Interface
266  Tag th;
267  CHKERR mwls_approx->mwlsMoab.tag_get_handle(
268  cp.mwlsStressTagName.c_str(), th);
269  CHKERR mwls_approx->getValuesToNodes(th);
270 
271  auto dm = simple_interface_ptr->getDM();
272  CHKERR DMoFEMLoopFiniteElements(
273  dm, simple_interface_ptr->getDomainFEName(), error_volume_fe);
274 
275  MWLSLog;
276  MOFEM_LOG_C("WORLD", Sev::inform, "Pre-process tag data < %s >",
277  name.c_str());
278  // get tag values at nodes
279  CHKERR mwls_approx->getValuesToNodes(tag_handles[t]);
280  // create tags on computational mesh for approximated values and
281  // derivatives of values
282  std::string approx_name = "APPROX_" + name;
283  std::string approx_diff_name[] = {"APPROX_DX_" + name,
284  "APPROX_DY_" + name,
285  "APPROX_DZ_" + name};
286  // get length and type of med tag on med mesh
287  int length;
288  CHKERR mwls_approx->mwlsMoab.tag_get_length(tag_handles[t],
289  length);
290  DataType data_type;
291  CHKERR mwls_approx->mwlsMoab.tag_get_data_type(tag_handles[t],
292  data_type);
293  // create tag on computational mesh
294  Tag th_approx;
295  CHKERR moab.tag_get_handle(
296  approx_name.c_str(), length, data_type, th_approx,
297  MB_TAG_CREAT | MB_TAG_SPARSE, def_vals);
298  // add tags to the lists
299  tag_approx_handles[t] = th_approx;
300  for (int d = 0; d != 3; d++) {
301  Tag th_approx_diff;
302  CHKERR moab.tag_get_handle(
303  approx_diff_name[d].c_str(), length, data_type,
304  th_approx_diff, MB_TAG_CREAT | MB_TAG_SPARSE, def_vals);
305  tag_approx_handles_diff(d, t) = th_approx_diff;
306  }
307 
308  std::string error_name = "STRESS_ERROR";
309  // create tag on computational mesh
310  Tag th_error;
311  CHKERR moab.tag_get_handle(error_name.c_str(), 1, data_type,
312  th_error, MB_TAG_CREAT | MB_TAG_SPARSE,
313  &def_error_vals);
314  // add tags to the lists
315  tag_error_handles[t] = th_error;
316 
317  std::string hydro_static_error_name = "HYDRO_STATIC_ERROR";
318  // create tag on computational mesh
319  Tag th_hydro_static_error;
320  CHKERR moab.tag_get_handle(hydro_static_error_name.c_str(), 1,
321  data_type, th_hydro_static_error,
322  MB_TAG_CREAT | MB_TAG_SPARSE,
323  &def_hydro_static_error_vals);
324  // add tags to the lists
325  tag_hydro_static_error_handles[t] = th_hydro_static_error;
326 
327  std::string deviatoric_error_name = "DEVIATORIC_ERROR";
328  // create tag on computational mesh
329  Tag th_deviatoric_error;
330  CHKERR moab.tag_get_handle(deviatoric_error_name.c_str(), 1,
331  data_type, th_deviatoric_error,
332  MB_TAG_CREAT | MB_TAG_SPARSE,
333  &def_deviatoric_error_vals);
334  // add tags to the lists
335  tag_deviatoric_error_handles[t] = th_deviatoric_error;
336  }
337  }
338 
339  // setup of the tags
340  if (name.compare("RHO") == 0) {
341  MWLSLog;
342  MOFEM_LOG_C("WORLD", Sev::inform, "Pre-process tag data < %s >",
343  name.c_str());
344 
345  CHKERR mwls_approx->getValuesToNodes(tag_handles[t]);
346  std::string approx_name = "RHO";
347  std::string approx_diff_name[] = {
348  "APPROX_DX_" + name, "APPROX_DY_" + name, "APPROX_DZ_" + name};
349  std::string approx_diff_diff_name[] = {
350  "APPROX_DX_DX_" + name, "APPROX_DX_DY_" + name,
351  "APPROX_DX_DZ_" + name, "APPROX_DY_DX_" + name,
352  "APPROX_DY_DY_" + name, "APPROX_DY_DZ_" + name,
353  "APPROX_DZ_DX_" + name, "APPROX_DZ_DY_" + name,
354  "APPROX_DZ_DZ_" + name};
355  int length;
356  CHKERR mwls_approx->mwlsMoab.tag_get_length(tag_handles[t], length);
357  DataType data_type;
358  CHKERR mwls_approx->mwlsMoab.tag_get_data_type(tag_handles[t],
359  data_type);
360  // create tag on computational mesh
361  Tag th_approx;
362  CHKERR moab.tag_get_handle(approx_name.c_str(), length, data_type,
363  th_approx, MB_TAG_CREAT | MB_TAG_SPARSE,
364  def_vals);
365  tag_approx_handles[t] = th_approx;
366  // add tags to the lists
367  for (int d = 0; d != 3; d++) {
368  Tag th_approx_diff;
369  CHKERR moab.tag_get_handle(
370  approx_diff_name[d].c_str(), length, data_type,
371  th_approx_diff, MB_TAG_CREAT | MB_TAG_SPARSE, def_vals);
372  tag_approx_handles_diff(d, t) = th_approx_diff;
373  for (int n = 0; n != 3; n++) {
374  Tag th_approx_diff_diff;
375  CHKERR moab.tag_get_handle(
376  approx_diff_diff_name[d + 3 * n].c_str(), length, data_type,
377  th_approx_diff_diff, MB_TAG_CREAT | MB_TAG_SPARSE,
378  def_vals);
379  tag_approx_handles_diff_diff(d + 3 * n, t) =
380  th_approx_diff_diff;
381  }
382  }
383  }
384  }
385 
386  EntityHandle part_meshset;
387  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, part_meshset);
388 
389  auto approx_and_eval_errors = [&](auto &tets) {
390  MoFEMFunctionBegin;
391 
392  double stress_sum = 0;
393  double stress_error_sum = 0;
394  double hydro_static_error_sum = 0;
395  double deviatoric_error_sum = 0;
396  double mesh_volume = 0;
397  double energy_sum = 0;
398  double energy_error_sum = 0;
399 
400  std::string printing_name;
401 
402  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
403 
404  auto mwls_approximate = [&](auto &tets) {
405  MoFEMFunctionBegin;
406 
407  auto fe_error_eval = boost::make_shared<FEMethod>();
408 
409  auto eval_error_op = [&]() -> MoFEMErrorCode {
410  MoFEMFunctionBegin;
411 
412  auto tet = fe_error_eval->numeredEntFiniteElementPtr->getEnt();
413  int rank = pcomm->rank();
414  CHKERR moab.tag_set_data(pcomm->part_tag(), &tet, 1, &rank);
415  CHKERR moab.add_entities(part_meshset, &tet, 1);
416  MOFEM_LOG("WORLD", Sev::noisy) << tet << " partition " << rank;
417 
418  int num_nodes;
419  const EntityHandle *conn;
420  CHKERR moab.get_connectivity(tet, conn, num_nodes, true);
421  MatrixDouble tet_coords(num_nodes, 3);
422  CHKERR moab.get_coords(conn, num_nodes,
423  &*tet_coords.data().begin());
424  const double vol = Tools::tetVolume(&*tet_coords.data().begin());
425  mesh_volume += vol;
426 
427  // get coords
428  double coords[3] = {0, 0, 0};
429  for (size_t n = 0; n != num_nodes; ++n)
430  for (auto d : {0, 1, 2})
431  coords[d] += tet_coords(n, d) / num_nodes;
432 
433  // find nodes in influence radius on med mesh
434  CHKERR mwls_approx->getInfluenceNodes(coords);
435  // calculate approximation/base functions
436  CHKERR mwls_approx->calculateApproxCoeff(true, true);
437  CHKERR mwls_approx->evaluateApproxFun(coords);
438  CHKERR mwls_approx->evaluateFirstDiffApproxFun(coords, true);
439  CHKERR mwls_approx->evaluateSecondDiffApproxFun(coords, false);
440 
441  // loop over all tags and use mwls to approximate med tags on
442  // computational nodes results save on the compositional mesh
443  // tags
444  MWLSLog;
445  for (int t = 0; t != tag_handles.size(); t++) {
446  std::string name;
447  CHKERR mwls_approx->mwlsMoab.tag_get_name(tag_handles[t], name);
448  if (name.compare(cp.mwlsStressTagName) == 0) {
449  printing_name = name;
450  if (tet == tets[0])
451  MOFEM_LOG_C("WORLD", Sev::inform, "Approximate tag < %s >",
452  name.c_str());
453 
454  // approximate data
455  CHKERR mwls_approx->getTagData(tag_handles[t]);
456  // get values
457  const auto &vals = mwls_approx->getDataApprox();
458  // get direvatives
459  const auto &diff_vals = mwls_approx->getDiffDataApprox();
460 
461  double total_energy_value = 0;
462  double total_energy_error_value = 0;
463  double total_stress_value = 0;
464  double total_stress_error_value = 0;
465  double hydro_static_error_value = 0;
466  double deviatoric_error_value = 0;
467 
468  CHKERR mwls_approx->getErrorAtNode(
469  tag_handles[t], total_stress_value,
470  total_stress_error_value, hydro_static_error_value,
471  deviatoric_error_value, total_energy_value,
472  total_energy_error_value, young_modulus, poisson_ratio);
473 
474  // save values on tags
475  CHKERR moab.tag_set_data(tag_approx_handles[t], &tet, 1,
476  &*vals.data().begin());
477  // save direvatives on tags
478  for (int d = 0; d != 3; d++)
479  CHKERR moab.tag_set_data(tag_approx_handles_diff(d, t),
480  &tet, 1, &diff_vals(d, 0));
481 
482  // save stress norm error
483  total_stress_error_value =
484  sqrt(total_stress_error_value / total_stress_value);
485  CHKERR moab.tag_set_data(tag_error_handles[t], &tet, 1,
486  &total_stress_error_value);
487 
488  // save hydro static stress on tags
489  hydro_static_error_value =
490  sqrt(hydro_static_error_value / total_stress_value);
491  CHKERR moab.tag_set_data(tag_hydro_static_error_handles[t],
492  &tet, 1, &hydro_static_error_value);
493 
494  // save deviatoric stress on tags
495  deviatoric_error_value =
496  sqrt(deviatoric_error_value / total_stress_value);
497  CHKERR moab.tag_set_data(tag_deviatoric_error_handles[t],
498  &tet, 1, &deviatoric_error_value);
499 
500  energy_sum += vol * total_energy_value;
501  energy_error_sum += vol * total_energy_error_value;
502  stress_sum += vol * total_stress_value;
503  stress_error_sum += vol * total_stress_error_value;
504  hydro_static_error_sum += vol * hydro_static_error_value;
505  deviatoric_error_sum += vol * deviatoric_error_value;
506  }
507 
508  if (name.compare("RHO") == 0) {
509 
510  if (tet == tets[0])
511  MOFEM_LOG_C("WORLD", Sev::inform, "Approximate tag < %s >",
512  name.c_str());
513 
514  // approximate data
515  CHKERR mwls_approx->getTagData(tag_handles[t]);
516 
517  // get values
518  const auto &vals = mwls_approx->getDataApprox();
519  // get direvatives
520  const auto &diff_vals = mwls_approx->getDiffDataApprox();
521  const auto &diff_diff_vals =
522  mwls_approx->getDiffDiffDataApprox();
523  // save values on tags
524  CHKERR moab.tag_set_data(tag_approx_handles[t], &tet, 1,
525  &*vals.data().begin());
526  // save direvatives on tags
527  for (int d = 0; d != 3; d++) {
528  CHKERR moab.tag_set_data(tag_approx_handles_diff(d, t),
529  &tet, 1, &diff_vals(d, 0));
530  for (int n = 0; n != 3; n++) {
531  CHKERR moab.tag_set_data(
532  tag_approx_handles_diff_diff(d + 3 * n, t), &tet, 1,
533  &diff_diff_vals(d + 3 * n, 0));
534  }
535  }
536  }
537  }
538 
539  MoFEMFunctionReturn(0);
540  };
541 
542  auto pre_proc = [&]() -> MoFEMErrorCode {
543  MoFEMFunctionBegin;
544  stress_sum = 0;
545  stress_error_sum = 0;
546  hydro_static_error_sum = 0;
547  deviatoric_error_sum = 0;
548  mesh_volume = 0;
549  energy_sum = 0;
550  energy_error_sum = 0;
551  MoFEMFunctionReturn(0);
552  };
553 
554  auto post_proc = [&]() -> MoFEMErrorCode {
555  MoFEMFunctionBegin;
556  auto petsc_vec = createSmartVectorMPI(
557  PETSC_COMM_WORLD, (!m_field.get_comm_rank()) ? 7 : 0, 7);
558  std::array<double, 7> data = {stress_sum,
559  stress_error_sum,
560  hydro_static_error_sum,
561  deviatoric_error_sum,
562  mesh_volume,
563  energy_sum,
564  energy_error_sum};
565  constexpr std::array<int, 7> indices = {0, 1, 2, 3, 4, 5, 6};
566  CHKERR VecSetValues(petsc_vec, 7, indices.data(), data.data(),
567  ADD_VALUES);
568  CHKERR VecAssemblyBegin(petsc_vec);
569  CHKERR VecAssemblyEnd(petsc_vec);
570  if (!m_field.get_comm_rank()) {
571  CHKERR VecGetValues(petsc_vec, 7, indices.data(), data.data());
572  stress_sum = data[0];
573  stress_error_sum = data[1];
574  hydro_static_error_sum = data[2];
575  deviatoric_error_sum = data[3];
576  mesh_volume = data[4];
577  energy_sum = data[5];
578  energy_error_sum = data[6];
579  }
580  MoFEMFunctionReturn(0);
581  };
582 
583  struct FEEvalError : public FEMethod {
584  FEEvalError() = default;
585  };
586 
587  fe_error_eval->preProcessHook = pre_proc;
588  fe_error_eval->operatorHook = eval_error_op;
589  fe_error_eval->postProcessHook = post_proc;
590 
591  auto dm = simple_interface_ptr->getDM();
592  CHKERR DMoFEMLoopFiniteElements(
593  dm, simple_interface_ptr->getDomainFEName(), fe_error_eval);
594 
595  stress_sum = sqrt(stress_sum);
596  stress_error_sum = sqrt(stress_error_sum);
597  hydro_static_error_sum = sqrt(hydro_static_error_sum);
598  deviatoric_error_sum = sqrt(deviatoric_error_sum);
599 
600  MoFEMFunctionReturn(0);
601  };
602 
603  auto print_information_out = [&]() {
604  MoFEMFunctionBegin;
605  MOFEM_LOG_TAG("WORLD", "mwls_at_gauss_point")
606  MOFEM_LOG_C("WORLD", Sev::inform, "Number of tetrahedrons < %d >",
607  tets.size());
608  if (printing_name.compare(cp.mwlsStressTagName) == 0) {
609  MOFEM_LOG_C("WORLD", Sev::inform, "Mesh volume < %e > ",
610  mesh_volume);
611  MOFEM_LOG_C("WORLD", Sev::inform, "Stress L2 norm < %e >",
612  stress_sum);
613  MOFEM_LOG_C(
614  "WORLD", Sev::inform,
615  "Stress error L2 norm <absolute> < %e > (relative) ( %e ) ",
616  stress_error_sum, stress_error_sum / stress_sum);
617  MOFEM_LOG_C("WORLD", Sev::inform,
618  "Hydrostatic stress error L2 norm <absolute> < %e "
619  "> (relative) ( %e ) ",
620  hydro_static_error_sum,
621  hydro_static_error_sum / stress_sum);
622  MOFEM_LOG_C(
623  "WORLD", Sev::inform,
624  "Deviator error L2 norm <absolute> < %e > (relative) ( %e ) ",
625  deviatoric_error_sum, deviatoric_error_sum / stress_sum);
626  MOFEM_LOG_C("WORLD", Sev::inform,
627  "Energy norm <absolute> < %e > ", energy_sum);
628  MOFEM_LOG_C(
629  "WORLD", Sev::inform,
630  "Energy norm error <absolute> < %e > (relative) ( %e ) ",
631  energy_error_sum, energy_error_sum / energy_sum);
632  }
633  MoFEMFunctionReturn(0);
634  };
635 
636  CHKERR mwls_approximate(tets);
637  CHKERR print_information_out();
638 
639  if (test == PETSC_TRUE) {
640  constexpr double test_tol = 1e-6;
641  if ((stress_error_sum / stress_sum) > test_tol)
642  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
643  "Approximation of stress has big error %e",
644  stress_error_sum / stress_sum);
645  }
646 
647  MoFEMFunctionReturn(0);
648  };
649 
650  // get nodes on tets on computational mesh
651  CHKERR approx_and_eval_errors(tets);
652 
653  // save meshes
654  CHKERR moab.write_file("out_mwls.h5m", "MOAB", "PARALLEL=WRITE_PART",
655  &part_meshset, 1);
656  }
657  }
658  }
659  CATCH_ERRORS;
660 
661  MoFEM::Core::Finalize();
662 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
FractureMechanics::CrackPropagation::residualStressBlock
int residualStressBlock
Block on which residual stress is applied.
Definition: CrackPropagation.hpp:136
EntityHandle
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
NOISY
@ NOISY
Definition: definitions.h:224
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MWLS.hpp
Mortar.hpp
MoFEM::Mat_Elastic
Elastic material data structure.
Definition: MaterialBlocks.hpp:139
MoFEM::createSmartVectorMPI
DEPRECATED auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Definition: PetscSmartObj.hpp:209
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM.hpp
MoFEM::BasicMethod::preProcessHook
boost::function< MoFEMErrorCode()> preProcessHook
Hook function for pre-processing.
Definition: LoopMethods.hpp:294
BasicFiniteElements.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::MeshsetsManager::getMeshset
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:708
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
FractureMechanics::MWLSApprox::OpMWLSStressAtGaussPts
Evaluate stress at integration points.
Definition: MWLS.hpp:404
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:747
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
FractureMechanics::ELASTIC_TAG
@ ELASTIC_TAG
Definition: CrackPropagation.hpp:37
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::MeshsetsManager::getEntitiesByDimension
MoFEMErrorCode getEntitiesByDimension(const int ms_id, const unsigned int cubit_bc_type, const int dimension, Range &entities, const bool recursive=true) const
get entities from CUBIT/meshset of a particular entity dimension
Definition: MeshsetsManager.cpp:669
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
Hooke.hpp
Implementation of linear elastic material.
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:172
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:368
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
CrackFrontElement.hpp
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:545
convert.n
n
Definition: convert.py:82
MoFEM::Tools::tetVolume
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition: Tools.cpp:30
MWLSLog
#define MWLSLog
Definition: mwls_at_gauss_point.cpp:45
Range
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
FractureMechanics::CrackPropagation::mwlsApproxFile
std::string mwlsApproxFile
Name of file with internal stresses.
Definition: CrackPropagation.hpp:132
main
int main(int argc, char *argv[])
Definition: mwls_at_gauss_point.cpp:60
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
FractureMechanics::CrackPropagation
Definition: CrackPropagation.hpp:77
FractureMechanics::CrackPropagation::mwlsStressTagName
std::string mwlsStressTagName
Name of tag with internal stresses.
Definition: CrackPropagation.hpp:133
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
help
static char help[]
Definition: mwls_at_gauss_point.cpp:42
ComplexConstArea.hpp
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::Simple::buildFields
MoFEMErrorCode buildFields()
Build fields.
Definition: Simple.cpp:554
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
FractureMechanics::MWLSApprox::OpMWLSStressAndErrorsAtGaussPts
Evaluate stress at integration points.
Definition: MWLS.hpp:454
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
CrackPropagation.hpp
Main class for crack propagation.
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
NeoHookean.hpp
Implementation of Neo-Hookean elastic material.
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:683
FractureMechanics::CrackFrontElement
CrackFrontSingularBase< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore > CrackFrontElement
Definition: CrackFrontElement.hpp:402
GriffithForceElement.hpp
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
FractureMechanics::CrackPropagation::getOptions
MoFEMErrorCode getOptions()
Get options form command line.
Definition: CrackPropagation.cpp:588
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
FractureMechanics
Definition: AnalyticalFun.hpp:15
Generated by Doxygen 1.8.17 and hosted at University of Glasgow