60 {
61 const char param_file[] = "param_file.petsc";
63
64 try {
65
67 PetscBool test;
70
71 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"mwls_at_gauss_point",
72 "none");
74 PETSC_NULL);
75 CHKERR PetscOptionsBool(
"-test",
"test code",
"", test, &test, PETSC_NULL);
76 CHKERR PetscOptionsScalar(
"-young_modulus",
77 "fix nodes below given threshold", "",
79 CHKERR PetscOptionsScalar(
"-poisson_ratio",
80 "fix nodes below given threshold", "",
82
83 ierr = PetscOptionsEnd();
85
86
87 DMType dm_name = "DMMOFEM";
89
90 moab::Core mb_instance;
91 moab::Interface &moab = mb_instance;
92
95
96
97 Simple *simple_interface_ptr;
99
101
102
107 bitLevel.set(0);
109 bitLevel);
110
111 CHKERR m_field.
get_moab().get_entities_by_dimension(whole_mesh, 3, ents,
112 true);
114 false);
115
116
121
125
127 "SPATIAL_POSITION");
129
131 "MESH_NODE_POSITIONS");
132
133 CHKERR m_field.
loop_dofs(
"MESH_NODE_POSITIONS", ent_method_material2, 0);
135
136
137 auto elastic_fe =
138 boost::make_shared<NonlinearElasticElement>(m_field,
ELASTIC_TAG);
139
140 {
141
144
146 CHKERR it->getAttributeDataStructure(mydata);
147 int id = it->getMeshsetId();
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
163
166 auto mwls_approx = boost::make_shared<MWLSApprox>(m_field);
167
168 {
169
170 CHKERR mwls_approx->loadMWLSMesh(cp.mwlsApproxFile.c_str());
171
175 if (cp.residualStressBlock != -1) {
177 approx_meshset);
178 }
179
181 CHKERR m_field.
get_moab().get_entities_by_type(whole_mesh, MBTET, tets,
182 true);
186
187
188 std::vector<Tag> tag_handles;
189
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
217 for (
int t = 0;
t != tag_handles.size();
t++) {
218
219 std::string name;
220 CHKERR mwls_approx->mwlsMoab.tag_get_name(tag_handles[
t], name);
221
222 if (name.compare(cp.mwlsStressTagName) == 0 &&
223 name.compare("RHO") != 0) {
224
225
227 Range ref_node_edges;
229 PetscBool boolean = PETSC_FALSE;
230
231 auto error_volume_fe = boost::shared_ptr<CrackFrontElement>(
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(
242 error_volume_fe->getOpPtrVector().push_back(
244 mat_pos_at_pts_ptr));
245
246 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr(
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(
257
258 error_volume_fe->getOpPtrVector().push_back(
260 mwls_approx, cp.mwlsStressTagName, m_field));
261
263 cp.residualStressBlock,
BLOCKSET, 3, mwls_approx->tetsInBlock,
264 true);
265
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();
274
276 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Pre-process tag data < %s >",
277 name.c_str());
278
279 CHKERR mwls_approx->getValuesToNodes(tag_handles[
t]);
280
281
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
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
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
299 tag_approx_handles[
t] = th_approx;
300 for (
int d = 0;
d != 3;
d++) {
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
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
315 tag_error_handles[
t] = th_error;
316
317 std::string hydro_static_error_name = "HYDRO_STATIC_ERROR";
318
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
325 tag_hydro_static_error_handles[
t] = th_hydro_static_error;
326
327 std::string deviatoric_error_name = "DEVIATORIC_ERROR";
328
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
335 tag_deviatoric_error_handles[
t] = th_deviatoric_error;
336 }
337 }
338
339
340 if (name.compare("RHO") == 0) {
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
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
367 for (
int d = 0;
d != 3;
d++) {
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
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
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;
420 CHKERR moab.get_connectivity(tet, conn, num_nodes,
true);
422 CHKERR moab.get_coords(conn, num_nodes,
423 &*tet_coords.data().begin());
425 mesh_volume += vol;
426
427
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
434 CHKERR mwls_approx->getInfluenceNodes(coords);
435
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
442
443
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
455 CHKERR mwls_approx->getTagData(tag_handles[
t]);
456
457 const auto &vals = mwls_approx->getDataApprox();
458
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,
473
474
475 CHKERR moab.tag_set_data(tag_approx_handles[
t], &tet, 1,
476 &*vals.data().begin());
477
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
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
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
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
515 CHKERR mwls_approx->getTagData(tag_handles[
t]);
516
517
518 const auto &vals = mwls_approx->getDataApprox();
519
520 const auto &diff_vals = mwls_approx->getDiffDataApprox();
521 const auto &diff_diff_vals =
522 mwls_approx->getDiffDiffDataApprox();
523
524 CHKERR moab.tag_set_data(tag_approx_handles[
t], &tet, 1,
525 &*vals.data().begin());
526
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++) {
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
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
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};
567 ADD_VALUES);
568 CHKERR VecAssemblyBegin(petsc_vec);
569 CHKERR VecAssemblyEnd(petsc_vec);
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();
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 = [&]() {
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);
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);
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)
643 "Approximation of stress has big error %e",
644 stress_error_sum / stress_sum);
645 }
646
648 };
649
650
651 CHKERR approx_and_eval_errors(tets);
652
653
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,...)
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#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"
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
#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
CrackFrontSingularBase< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore > CrackFrontElement
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
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.
double poisson_ratio
Poisson ratio.
constexpr double t
plate stiffness
Evaluate stress at integration points.
Evaluate stress at integration points.
virtual moab::Interface & get_moab()=0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
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.
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.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode buildFields()
Build fields.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getDomainFEName() const
Get the Domain FE Name.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.