v0.14.0
Loading...
Searching...
No Matches
mwls_at_gauss_point.cpp File Reference
#include <tetgen.h>
#include <moab/SpatialLocator.hpp>
#include <MoFEM.hpp>
#include <BasicFiniteElements.hpp>
#include <Mortar.hpp>
#include <NeoHookean.hpp>
#include <Hooke.hpp>
#include <CrackFrontElement.hpp>
#include <ComplexConstArea.hpp>
#include <MWLS.hpp>
#include <GriffithForceElement.hpp>
#include <VolumeLengthQuality.hpp>
#include <CrackPropagation.hpp>

Go to the source code of this file.

Macros

#define BOOST_UBLAS_NDEBUG   1
 
#define MWLSLog
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help []
 

Macro Definition Documentation

◆ BOOST_UBLAS_NDEBUG

#define BOOST_UBLAS_NDEBUG   1

Definition at line 32 of file mwls_at_gauss_point.cpp.

◆ MWLSLog

#define MWLSLog
Value:
MOFEM_LOG_CHANNEL("WORLD"); \
MOFEM_LOG_FUNCTION(); \
MOFEM_LOG_TAG("WORLD", "mwls_at_gauss_point");
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.

Definition at line 45 of file mwls_at_gauss_point.cpp.

45#define MWLSLog \
46 MOFEM_LOG_CHANNEL("WORLD"); \
47 MOFEM_LOG_FUNCTION(); \
48 MOFEM_LOG_TAG("WORLD", "mwls_at_gauss_point");

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 60 of file mwls_at_gauss_point.cpp.

60 {
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();
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,
119 CHKERR simple_interface_ptr->addDomainField("SPATIAL_POSITION", H1,
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
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(
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(
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();
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) {
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) {
406
407 auto fe_error_eval = boost::make_shared<FEMethod>();
408
409 auto eval_error_op = [&]() -> MoFEMErrorCode {
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
540 };
541
542 auto pre_proc = [&]() -> MoFEMErrorCode {
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;
552 };
553
554 auto post_proc = [&]() -> MoFEMErrorCode {
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 }
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();
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
601 };
602
603 auto print_information_out = [&]() {
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);
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);
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);
629 "WORLD", Sev::inform,
630 "Energy norm error <absolute> < %e > (relative) ( %e ) ",
631 energy_error_sum, energy_error_sum / energy_sum);
632 }
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
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 }
660
662}
#define MOFEM_LOG_C(channel, severity, format,...)
@ NOISY
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ H1
continuous field
Definition definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
@ BLOCKSET
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
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
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
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.
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
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
const double n
refractive index of diffusive medium
static char help[]
#define MWLSLog
CrackFrontSingularBase< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore > CrackFrontElement
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
UBlasVector< double > VectorDouble
Definition Types.hpp:68
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
DEPRECATED auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
double young_modulus
Young modulus.
Definition plastic.cpp:121
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:122
constexpr double t
plate stiffness
Definition plate.cpp:58
Evaluate stress at integration points.
Definition MWLS.hpp:405
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
Deprecated interface functions.
structure for User Loop Methods on finite elements
Elastic material data structure.
Interface for managing meshsets containing materials and boundary conditions.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
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
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:194
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:829
MoFEMErrorCode buildFields()
Build fields.
Definition Simple.cpp:608
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:599
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:765
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:387
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition Tools.cpp:30
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.

Variable Documentation

◆ help

char help[]
static
Initial value:
= "Testing moving weighted least square approximation"
"n\n"

Definition at line 42 of file mwls_at_gauss_point.cpp.