v0.13.1
MixTransportElement.hpp
Go to the documentation of this file.
1/** \file MixTransportElement.hpp
2 * \brief Mix implementation of transport element
3 *
4 * \ingroup mofem_mix_transport_elem
5 *
6 */
7
8
9
10#ifndef _MIX_TRANPORT_ELEMENT_HPP_
11#define _MIX_TRANPORT_ELEMENT_HPP_
12
13namespace MixTransport {
14
15/** \brief Mix transport problem
16 \ingroup mofem_mix_transport_elem
17
18 Note to solve this system you need to use direct solver or proper
19 preconditioner for saddle problem.
20
21 It is based on \cite arnold2006differential \cite arnold2012mixed \cite
22 reviartthomas1996
23 <https://www.researchgate.net/profile/Richard_Falk/publication/226454406_Differential_Complexes_and_Stability_of_Finite_Element_Methods_I._The_de_Rham_Complex/links/02e7e5214f0426ff77000000.pdf>
24
25 General problem have form,
26 \f[
27 \mathbf{A} \boldsymbol\sigma + \textrm{grad}[u] = \mathbf{0} \; \textrm{on} \;
28 \Omega \\ \textrm{div}[\boldsymbol\sigma] = f \; \textrm{on} \; \Omega \f]
29
30*/
32
34
35 /**
36 * \brief definition of volume element
37
38 * It is used to calculate volume integrals. On volume element we set-up
39 * operators to calculate components of matrix and vector.
40
41 */
45 int getRule(int order) { return 2 * order + 1; };
46 };
47
48 MyVolumeFE feVol; ///< Instance of volume element
49
50 /** \brief definition of surface element
51
52 * It is used to calculate surface integrals. On volume element are operators
53 * evaluating natural boundary conditions.
54
55 */
59 int getRule(int order) { return 2 * order + 1; };
60 };
61
62 MyTriFE feTri; ///< Instance of surface element
63
64 /**
65 * \brief construction of this data structure
66 */
68 : mField(m_field), feVol(m_field), feTri(m_field){};
69
70 /**
71 * \brief destructor
72 */
74
75 VectorDouble valuesAtGaussPts; ///< values at integration points on element
77 valuesGradientAtGaussPts; ///< gradients at integration points on element
79 divergenceAtGaussPts; ///< divergence at integration points on element
80 MatrixDouble fluxesAtGaussPts; ///< fluxes at integration points on element
81
82 set<int> bcIndices;
83
84 /**
85 * \brief get dof indices where essential boundary conditions are applied
86 * @param is indices
87 * @return error code
88 */
91 std::vector<int> ids;
92 ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
93 IS is_local;
94 CHKERR ISCreateGeneral(mField.get_comm(), ids.size(),
95 ids.empty() ? PETSC_NULL : &ids[0],
96 PETSC_COPY_VALUES, &is_local);
97 CHKERR ISAllGather(is_local, is);
98 CHKERR ISDestroy(&is_local);
100 }
101
102 /**
103 * \brief set source term
104 * @param ent handle to entity on which function is evaluated
105 * @param x coord
106 * @param y coord
107 * @param z coord
108 * @param flux reference to source term set by function
109 * @return error code
110 */
111 virtual MoFEMErrorCode getSource(const EntityHandle ent, const double x,
112 const double y, const double z,
113 double &flux) {
115 flux = 0;
117 }
118
119 /**
120 * \brief natural (Dirichlet) boundary conditions (set values)
121 * @param ent handle to finite element entity
122 * @param x coord
123 * @param y coord
124 * @param z coord
125 * @param value reference to value set by function
126 * @return error code
127 */
128 virtual MoFEMErrorCode getResistivity(const EntityHandle ent, const double x,
129 const double y, const double z,
130 MatrixDouble3by3 &inv_k) {
132 inv_k.clear();
133 for (int dd = 0; dd < 3; dd++) {
134 inv_k(dd, dd) = 1;
135 }
137 }
138
139 /**
140 * \brief evaluate natural (Dirichlet) boundary conditions
141 * @param ent entity on which bc is evaluated
142 * @param x coordinate
143 * @param y coordinate
144 * @param z coordinate
145 * @param value vale
146 * @return error code
147 */
148 virtual MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg,
149 const double x, const double y,
150 const double z, double &value) {
152 value = 0;
154 }
155
156 /**
157 * \brief essential (Neumann) boundary condition (set fluxes)
158 * @param ent handle to finite element entity
159 * @param x coord
160 * @param y coord
161 * @param z coord
162 * @param flux reference to flux which is set by function
163 * @return [description]
164 */
165 virtual MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x,
166 const double y, const double z,
167 double &flux) {
169 flux = 0;
171 }
172
173 /** \brief data for evaluation of het conductivity and heat capacity elements
174 * \ingroup mofem_thermal_elem
175 */
176 struct BlockData {
178 double cApacity;
179 Range tEts; ///< constatins elements in block set
180 };
181 std::map<int, BlockData>
182 setOfBlocks; ///< maps block set id with appropriate BlockData
183
184 /**
185 * \brief Add fields to database
186 * @param values name of the fields
187 * @param fluxes name of filed for fluxes
188 * @param order order of approximation
189 * @return error code
190 */
191 MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes,
192 const int order) {
194 // Fields
197
198 // meshset consisting all entities in mesh
199 EntityHandle root_set = mField.get_moab().get_root_set();
200 // add entities to field
201 CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, fluxes);
202 CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, values);
203 CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
204 CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
205 CHKERR mField.set_field_order(root_set, MBTET, values, order);
207 }
208
209 /// \brief add finite elements
211 const std::string &fluxes_name, const std::string &values_name) {
213
214 // Set up volume element operators. Operators are used to calculate
215 // components of stiffness matrix & right hand side, in essence are used to
216 // do volume integrals over tetrahedral in this case.
217
218 // Define element "MIX". Note that this element will work with fluxes_name
219 // and values_name. This reflect bilinear form for the problem
227
228 // Define finite element to integrate over skeleton, we need that to
229 // evaluate error
230 CHKERR mField.add_finite_element("MIX_SKELETON", MF_ZERO);
232 fluxes_name);
234 fluxes_name);
236 fluxes_name);
237
238 // Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those
239 // BLOCKSETS and add them to "MIX" finite element. In addition get data form
240 // that meshset and set conductivity which is used to calculate fluxes from
241 // gradients of concentration or gradient of temperature, depending how you
242 // interpret variables.
245
246 Mat_Thermal temp_data;
247 CHKERR it->getAttributeDataStructure(temp_data);
248 setOfBlocks[it->getMeshsetId()].cOnductivity =
249 temp_data.data.Conductivity;
250 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
251 CHKERR mField.get_moab().get_entities_by_type(
252 it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
254 setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
255
256 Range skeleton;
257 CHKERR mField.get_moab().get_adjacencies(
258 setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
259 moab::Interface::UNION);
261 "MIX_SKELETON");
262 }
263
264 // Define element to integrate natural boundary conditions, i.e. set values.
265 CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
267 fluxes_name);
269 fluxes_name);
271 fluxes_name);
273 values_name);
274
275 // Define element to apply essential boundary conditions.
278 fluxes_name);
280 fluxes_name);
282 fluxes_name);
284 values_name);
286 }
287
288 /**
289 * \brief Build problem
290 * @param ref_level mesh refinement on which mesh problem you like to built.
291 * @return error code
292 */
295 // build field
297 // get tetrahedrons which has been build previously and now in so called
298 // garbage bit level
299 Range done_tets;
300 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
301 BitRefLevel().set(0), BitRefLevel().set(), MBTET, done_tets);
302 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
303 BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTET,
304 done_tets);
305 // get tetrahedrons which belong to problem bit level
306 Range ref_tets;
307 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
308 ref_level, BitRefLevel().set(), MBTET, ref_tets);
309 ref_tets = subtract(ref_tets, done_tets);
310 CHKERR mField.build_finite_elements("MIX", &ref_tets, 2);
311 // get triangles which has been build previously and now in so called
312 // garbage bit level
313 Range done_faces;
314 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
315 BitRefLevel().set(0), BitRefLevel().set(), MBTRI, done_faces);
316 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
317 BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTRI,
318 done_faces);
319 // get triangles which belong to problem bit level
320 Range ref_faces;
321 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
322 ref_level, BitRefLevel().set(), MBTRI, ref_faces);
323 ref_faces = subtract(ref_faces, done_faces);
324 // build finite elements structures
325 CHKERR mField.build_finite_elements("MIX_BCFLUX", &ref_faces, 2);
326 CHKERR mField.build_finite_elements("MIX_BCVALUE", &ref_faces, 2);
327 CHKERR mField.build_finite_elements("MIX_SKELETON", &ref_faces, 2);
328 // Build adjacencies of degrees of freedom and elements
330 // Define problem
332 // set refinement level for problem
334 // Add element to problem
336 CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_SKELETON");
337 // Boundary conditions
339 CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCVALUE");
340 // build problem
341
342 ProblemsManager *prb_mng_ptr;
343 CHKERR mField.getInterface(prb_mng_ptr);
344 CHKERR prb_mng_ptr->buildProblem("MIX", true);
345 // mesh partitioning
346 // partition
347 CHKERR prb_mng_ptr->partitionProblem("MIX");
348 CHKERR prb_mng_ptr->partitionFiniteElements("MIX");
349 // what are ghost nodes, see Petsc Manual
350 CHKERR prb_mng_ptr->partitionGhostDofs("MIX");
352 }
353
357 std::vector<EntityHandle> &mapGaussPts;
359 std::vector<EntityHandle> &map_gauss_pts)
360 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
361 "VALUES", UserDataOperator::OPCOL),
362 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
366 if (type != MBTET)
368 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
369 Tag th_error_flux;
370 CHKERR getVolumeFE()->mField.get_moab().tag_get_handle("ERROR_FLUX",
371 th_error_flux);
372 double *error_flux_ptr;
373 CHKERR getVolumeFE()->mField.get_moab().tag_get_by_ptr(
374 th_error_flux, &fe_ent, 1, (const void **)&error_flux_ptr);
375
376 Tag th_error_div;
377 CHKERR getVolumeFE()->mField.get_moab().tag_get_handle("ERROR_DIV",
378 th_error_div);
379 double *error_div_ptr;
380 CHKERR getVolumeFE()->mField.get_moab().tag_get_by_ptr(
381 th_error_div, &fe_ent, 1, (const void **)&error_div_ptr);
382
383 Tag th_error_jump;
384 CHKERR getVolumeFE()->mField.get_moab().tag_get_handle("ERROR_JUMP",
385 th_error_jump);
386 double *error_jump_ptr;
387 CHKERR getVolumeFE()->mField.get_moab().tag_get_by_ptr(
388 th_error_jump, &fe_ent, 1, (const void **)&error_jump_ptr);
389
390 {
391 double def_val = 0;
392 Tag th_error_flux;
393 CHKERR postProcMesh.tag_get_handle(
394 "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
395 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
396 for (vector<EntityHandle>::iterator vit = mapGaussPts.begin();
397 vit != mapGaussPts.end(); vit++) {
398 CHKERR postProcMesh.tag_set_data(th_error_flux, &*vit, 1,
399 error_flux_ptr);
400 }
401
402 Tag th_error_div;
403 CHKERR postProcMesh.tag_get_handle(
404 "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
405 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
406 for (vector<EntityHandle>::iterator vit = mapGaussPts.begin();
407 vit != mapGaussPts.end(); vit++) {
408 CHKERR postProcMesh.tag_set_data(th_error_div, &*vit, 1,
409 error_div_ptr);
410 }
411
412 Tag th_error_jump;
413 CHKERR postProcMesh.tag_get_handle(
414 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
415 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
416 for (vector<EntityHandle>::iterator vit = mapGaussPts.begin();
417 vit != mapGaussPts.end(); vit++) {
418 CHKERR postProcMesh.tag_set_data(th_error_jump, &*vit, 1,
419 error_jump_ptr);
420 }
421 }
423 }
424 };
425
426 /**
427 * \brief Post process results
428 * @return error code
429 */
430 MoFEMErrorCode postProc(const string out_file) {
432
433 PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
434 mField);
435
436 auto jac_ptr = boost::make_shared<MatrixDouble>();
437 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
438 auto det_ptr = boost::make_shared<VectorDouble>();
439 post_proc.getOpPtrVector().push_back(new OpCalculateHOJac<3>(jac_ptr));
440 post_proc.getOpPtrVector().push_back(
441 new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
442 post_proc.getOpPtrVector().push_back(
443 new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
444 post_proc.getOpPtrVector().push_back(
445 new OpSetHOInvJacToScalarBases<3>(L2, inv_jac_ptr));
446
447 auto values_ptr = boost::make_shared<VectorDouble>();
448 auto grad_ptr = boost::make_shared<MatrixDouble>();
449 auto flux_ptr = boost::make_shared<MatrixDouble>();
450
451 post_proc.getOpPtrVector().push_back(
452 new OpCalculateScalarFieldValues("VALUES", values_ptr));
453 post_proc.getOpPtrVector().push_back(
454 new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
455 post_proc.getOpPtrVector().push_back(
456 new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
457
458 using OpPPMap = OpPostProcMapInMoab<3, 3>;
459
460 post_proc.getOpPtrVector().push_back(
461
462 new OpPPMap(post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
463
464 {{"VALUES", values_ptr}},
465
466 {{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
467
468 {}, {}
469
470 ));
471
472 post_proc.getOpPtrVector().push_back(new OpPostProc(
473 post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
474
475 CHKERR mField.loop_finite_elements("MIX", "MIX", post_proc);
476
477 CHKERR post_proc.writeFile(out_file.c_str());
479 }
480
481 Vec D, D0, F;
482 Mat Aij;
483
484 /// \brief create matrices
487 CHKERR mField.getInterface<MatrixManager>()
488 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("MIX", &Aij);
489 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
490 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
491 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
493 }
494
495 /**
496 * \brief solve problem
497 * @return error code
498 */
501 CHKERR MatZeroEntries(Aij);
502 CHKERR VecZeroEntries(F);
503 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
504 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
505 CHKERR VecZeroEntries(D0);
506 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
507 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
508 CHKERR VecZeroEntries(D);
509 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
510 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
511
512 CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
513 "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
514
515 // Calculate essential boundary conditions
516
517 // clear essential bc indices, it could have dofs from other mesh refinement
518 bcIndices.clear();
519 // clear operator, just in case if some other operators are left on this
520 // element
521 feTri.getOpPtrVector().clear();
522 feTri.getOpPtrVector().push_back(
523 new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
524 // set operator to calculate essential boundary conditions
525 feTri.getOpPtrVector().push_back(new OpEvaluateBcOnFluxes(*this, "FLUXES"));
526 CHKERR mField.loop_finite_elements("MIX", "MIX_BCFLUX", feTri);
527 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
528 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
529 CHKERR VecAssemblyBegin(D0);
530 CHKERR VecAssemblyEnd(D0);
531
532 // set operators to calculate matrix and right hand side vectors
533 feVol.getOpPtrVector().clear();
534 feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
535 feVol.getOpPtrVector().push_back(
536 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
537 feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
538 feVol.getOpPtrVector().push_back(
539 new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
540 feVol.getOpPtrVector().push_back(
541 new OpTauDotSigma_HdivHdiv(*this, "FLUXES", Aij, F));
542 feVol.getOpPtrVector().push_back(
543 new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", Aij, F));
544 CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
545
546 // calculate right hand side for natural boundary conditions
547 feTri.getOpPtrVector().clear();
548 feTri.getOpPtrVector().push_back(
549 new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
550 feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
551 CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
552
553 // assemble matrices
554 CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
555 CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
556 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
557 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
558 CHKERR VecAssemblyBegin(F);
559 CHKERR VecAssemblyEnd(F);
560
561 {
562 double nrm2_F;
563 CHKERR VecNorm(F, NORM_2, &nrm2_F);
564 PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
565 }
566
567 // CHKERR MatMultAdd(Aij,D0,F,F); ;
568
569 // for ksp solver vector is moved into rhs side
570 // for snes it is left ond the left
571 CHKERR VecScale(F, -1);
572
573 IS essential_bc_ids;
574 CHKERR getDirichletBCIndices(&essential_bc_ids);
575 CHKERR MatZeroRowsColumnsIS(Aij, essential_bc_ids, 1, D0, F);
576 CHKERR ISDestroy(&essential_bc_ids);
577
578 // {
579 // double norm;
580 // CHKERR MatNorm(Aij,NORM_FROBENIUS,&norm); ;
581 // PetscPrintf(PETSC_COMM_WORLD,"mat norm = %6.4e\n",norm);
582 // }
583
584 {
585 double nrm2_F;
586 CHKERR VecNorm(F, NORM_2, &nrm2_F);
587 PetscPrintf(PETSC_COMM_WORLD, "With BC nrm2_F = %6.4e\n", nrm2_F);
588 }
589
590 // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
591 // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
592 // std::string wait;
593 // std::cin >> wait;
594
595 // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
596 // std::cin >> wait;
597
598 // Solve
599 KSP solver;
600 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
601 CHKERR KSPSetOperators(solver, Aij, Aij);
602 CHKERR KSPSetFromOptions(solver);
603 CHKERR KSPSetUp(solver);
604 CHKERR KSPSolve(solver, F, D);
605 CHKERR KSPDestroy(&solver);
606
607 {
608 double nrm2_D;
609 CHKERR VecNorm(D, NORM_2, &nrm2_D);
610 ;
611 PetscPrintf(PETSC_COMM_WORLD, "nrm2_D = %6.4e\n", nrm2_D);
612 }
613 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
614 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
615
616 // copy data form vector on mesh
617 CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
618 "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
619
621 }
622
623 /// \brief calculate residual
626 CHKERR VecZeroEntries(F);
627 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
628 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
629 CHKERR VecAssemblyBegin(F);
630 CHKERR VecAssemblyEnd(F);
631 // calculate residuals
632 feVol.getOpPtrVector().clear();
633 feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
634 feVol.getOpPtrVector().push_back(
635 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
636 feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
637 feVol.getOpPtrVector().push_back(
638 new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
639 feVol.getOpPtrVector().push_back(
640 new OpTauDotSigma_HdivHdiv(*this, "FLUXES", PETSC_NULL, F));
641 feVol.getOpPtrVector().push_back(
642 new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", PETSC_NULL, F));
643 CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
644 feTri.getOpPtrVector().clear();
645 feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
646 CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
647 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
648 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
649 CHKERR VecAssemblyBegin(F);
650 CHKERR VecAssemblyEnd(F);
651 // CHKERR VecAXPY(F,-1.,D0);
652 // CHKERR MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULL,F);
653 {
654 std::vector<int> ids;
655 ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
656 std::vector<double> vals(ids.size(), 0);
657 CHKERR VecSetValues(F, ids.size(), &*ids.begin(), &*vals.begin(),
658 INSERT_VALUES);
659 CHKERR VecAssemblyBegin(F);
660 CHKERR VecAssemblyEnd(F);
661 }
662 {
663 double nrm2_F;
664 CHKERR VecNorm(F, NORM_2, &nrm2_F);
665 PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
666 const double eps = 1e-8;
667 if (nrm2_F > eps) {
668 // SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"problem with
669 // residual");
670 }
671 }
673 }
674
675 /** \brief Calculate error on elements
676
677 \todo this functions runs serial, in future should be parallel and work on
678 distributed meshes.
679
680 */
683 errorMap.clear();
684 sumErrorFlux = 0;
685 sumErrorDiv = 0;
686 sumErrorJump = 0;
687 feTri.getOpPtrVector().clear();
688 feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
689 CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
690 mField.get_comm_size());
691 feVol.getOpPtrVector().clear();
692 feVol.getOpPtrVector().push_back(
693 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
694 feVol.getOpPtrVector().push_back(
695 new OpValuesGradientAtGaussPts(*this, "VALUES"));
696 feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
697 CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
698 mField.get_comm_size());
699 const Problem *problem_ptr;
700 CHKERR mField.get_problem("MIX", &problem_ptr);
701 PetscPrintf(mField.get_comm(),
702 "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
703 "jump^2 = %6.4e error tot^2 = %6.4e\n",
704 problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
705 sumErrorJump, sumErrorFlux + sumErrorDiv + sumErrorJump);
707 }
708
709 /// \brief destroy matrices
712 CHKERR MatDestroy(&Aij);
713 ;
714 CHKERR VecDestroy(&D);
715 ;
716 CHKERR VecDestroy(&D0);
717 ;
718 CHKERR VecDestroy(&F);
719 ;
721 }
722
723 /**
724 \brief Assemble \f$\int_\mathcal{T} \mathbf{A} \boldsymbol\sigma \cdot
725 \boldsymbol\tau \textrm{d}\mathcal{T}\f$
726
727 \ingroup mofem_mix_transport_elem
728 */
731
733 Mat Aij;
735
737 const std::string flux_name, Mat aij, Vec f)
738 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
739 flux_name, flux_name,
740 UserDataOperator::OPROWCOL | UserDataOperator::OPCOL),
741 cTx(ctx), Aij(aij), F(f) {
742 sYmm = true;
743 }
744
745 MatrixDouble NN, transNN;
748
749 /**
750 * \brief Assemble matrix
751 * @param row_side local index of row entity on element
752 * @param col_side local index of col entity on element
753 * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
754 * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
755 * @param row_data data for row
756 * @param col_data data for col
757 * @return error code
758 */
759 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
760 EntityType col_type,
762 EntitiesFieldData::EntData &col_data) {
764 if (Aij == PETSC_NULL)
766 if (row_data.getIndices().size() == 0)
768 if (col_data.getIndices().size() == 0)
770 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
771 int nb_row = row_data.getIndices().size();
772 int nb_col = col_data.getIndices().size();
773 NN.resize(nb_row, nb_col, false);
774 NN.clear();
777 invK.resize(3, 3, false);
778 // get access to resistivity data by tensor rank 2
780 &invK(0, 0), &invK(0, 1), &invK(0, 2), &invK(1, 0), &invK(1, 1),
781 &invK(1, 2), &invK(2, 0), &invK(2, 1), &invK(2, 2));
782 // Get base functions
783 auto t_n_hdiv_row = row_data.getFTensor1N<3>();
785 int nb_gauss_pts = row_data.getN().size1();
786 for (int gg = 0; gg != nb_gauss_pts; gg++) {
787 // get integration weight and multiply by element volume
788 double w = getGaussPts()(3, gg) * getVolume();
789 const double x = getCoordsAtGaussPts()(gg, 0);
790 const double y = getCoordsAtGaussPts()(gg, 1);
791 const double z = getCoordsAtGaussPts()(gg, 2);
792 // calculate receptivity (invers of conductivity)
793 CHKERR cTx.getResistivity(fe_ent, x, y, z, invK);
794 for (int kk = 0; kk != nb_row; kk++) {
796 &col_data.getVectorN<3>(gg)(0, HVEC0),
797 &col_data.getVectorN<3>(gg)(0, HVEC1),
798 &col_data.getVectorN<3>(gg)(0, HVEC2), 3);
799 t_row(j) = w * t_n_hdiv_row(i) * t_inv_k(i, j);
800 for (int ll = 0; ll != nb_col; ll++) {
801 NN(kk, ll) += t_row(j) * t_n_hdiv_col(j);
802 ++t_n_hdiv_col;
803 }
804 ++t_n_hdiv_row;
805 }
806 }
807 Mat a = (Aij != PETSC_NULL) ? Aij : getFEMethod()->ts_B;
808 CHKERR MatSetValues(a, nb_row, &row_data.getIndices()[0], nb_col,
809 &col_data.getIndices()[0], &NN(0, 0), ADD_VALUES);
810 // matrix is symmetric, assemble other part
811 if (row_side != col_side || row_type != col_type) {
812 transNN.resize(nb_col, nb_row);
813 noalias(transNN) = trans(NN);
814 CHKERR MatSetValues(a, nb_col, &col_data.getIndices()[0], nb_row,
815 &row_data.getIndices()[0], &transNN(0, 0),
816 ADD_VALUES);
817 }
818
820 }
821
822 /**
823 * \brief Assemble matrix
824 * @param side local index of row entity on element
825 * @param type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
826 * @param data data for row
827 * @return error code
828 */
832 if (F == PETSC_NULL)
834 int nb_row = data.getIndices().size();
835 if (nb_row == 0)
837
838 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
839 Nf.resize(nb_row);
840 Nf.clear();
843 invK.resize(3, 3, false);
844 Nf.resize(nb_row);
845 Nf.clear();
846 // get access to resistivity data by tensor rank 2
848 &invK(0, 0), &invK(0, 1), &invK(0, 2), &invK(1, 0), &invK(1, 1),
849 &invK(1, 2), &invK(2, 0), &invK(2, 1), &invK(2, 2));
850 // get base functions
851 auto t_n_hdiv = data.getFTensor1N<3>();
852 auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
853 int nb_gauss_pts = data.getN().size1();
854 for (int gg = 0; gg < nb_gauss_pts; gg++) {
855 double w = getGaussPts()(3, gg) * getVolume();
856 const double x = getCoordsAtGaussPts()(gg, 0);
857 const double y = getCoordsAtGaussPts()(gg, 1);
858 const double z = getCoordsAtGaussPts()(gg, 2);
859 CHKERR cTx.getResistivity(fe_ent, x, y, z, invK);
860 for (int ll = 0; ll != nb_row; ll++) {
861 Nf[ll] += w * t_n_hdiv(i) * t_inv_k(i, j) * t_flux(j);
862 ++t_n_hdiv;
863 }
864 ++t_flux;
865 }
866 CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
867
869 }
870 };
871
872 /** \brief Assemble \f$ \int_\mathcal{T} u \textrm{div}[\boldsymbol\tau]
873 * \textrm{d}\mathcal{T} \f$
874 */
877
880
881 OpDivTauU_HdivL2(MixTransportElement &ctx, const std::string flux_name_row,
882 string val_name_col, Vec f)
883 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
884 flux_name_row, val_name_col, UserDataOperator::OPROW),
885 cTx(ctx), F(f) {
886 // this operator is not symmetric setting this variable makes element
887 // operator to loop over element entities (sub-simplices) without
888 // assumption that off-diagonal matrices are symmetric.
889 sYmm = false;
890 }
891
893
897 if (data.getFieldData().size() == 0)
899 int nb_row = data.getIndices().size();
900 Nf.resize(nb_row);
901 Nf.clear();
902 divVec.resize(data.getN().size2() / 3, 0);
903 if (divVec.size() != data.getIndices().size()) {
904 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
905 "data inconsistency");
906 }
907 int nb_gauss_pts = data.getN().size1();
908
910 auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
911
912 int gg = 0;
913 for (; gg < nb_gauss_pts; gg++) {
914 double w = getGaussPts()(3, gg) * getVolume();
915 for(auto &v : divVec) {
916 v = t_base_diff_hdiv(i, i);
917 ++t_base_diff_hdiv;
918 }
919 noalias(Nf) -= w * divVec * cTx.valuesAtGaussPts[gg];
920 }
921 CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
922 ;
923
925 }
926 };
927
928 /** \brief \f$ \int_\mathcal{T} \textrm{div}[\boldsymbol\sigma] v
929 * \textrm{d}\mathcal{T} \f$
930 */
933
935 Mat Aij;
937
938 /**
939 * \brief Constructor
940 */
941 OpVDivSigma_L2Hdiv(MixTransportElement &ctx, const std::string val_name_row,
942 string flux_name_col, Mat aij, Vec f)
943 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
944 val_name_row, flux_name_col,
945 UserDataOperator::OPROW | UserDataOperator::OPROWCOL),
946 cTx(ctx), Aij(aij), F(f) {
947
948 // this operator is not symmetric setting this variable makes element
949 // operator to loop over element entities without
950 // assumption that off-diagonal matrices are symmetric.
951 sYmm = false;
952 }
953
954 MatrixDouble NN, transNN;
956
957 /**
958 * \brief Do calculations
959 * @param row_side local index of entity on row
960 * @param col_side local index of entity on column
961 * @param row_type type of row entity
962 * @param col_type type of col entity
963 * @param row_data row data structure carrying information about base
964 * functions, DOFs indices, etc.
965 * @param col_data column data structure carrying information about base
966 * functions, DOFs indices, etc.
967 * @return error code
968 */
969 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
970 EntityType col_type,
972 EntitiesFieldData::EntData &col_data) {
974 if (Aij == PETSC_NULL)
976 if (row_data.getFieldData().size() == 0)
978 if (col_data.getFieldData().size() == 0)
980 int nb_row = row_data.getFieldData().size();
981 int nb_col = col_data.getFieldData().size();
982 NN.resize(nb_row, nb_col);
983 NN.clear();
984 divVec.resize(nb_col, false);
985
987 auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
988
989 int nb_gauss_pts = row_data.getN().size1();
990 for (int gg = 0; gg < nb_gauss_pts; gg++) {
991 double w = getGaussPts()(3, gg) * getVolume();
992 for (auto &v : divVec) {
993 v = t_base_diff_hdiv(i, i);
994 ++t_base_diff_hdiv;
995 }
996 noalias(NN) += w * outer_prod(row_data.getN(gg), divVec);
997 }
998 CHKERR MatSetValues(Aij, nb_row, &row_data.getIndices()[0], nb_col,
999 &col_data.getIndices()[0], &NN(0, 0), ADD_VALUES);
1000 transNN.resize(nb_col, nb_row);
1001 ublas::noalias(transNN) = -trans(NN);
1002 CHKERR MatSetValues(Aij, nb_col, &col_data.getIndices()[0], nb_row,
1003 &row_data.getIndices()[0], &transNN(0, 0),
1004 ADD_VALUES);
1006 }
1007
1011 if (data.getIndices().size() == 0)
1013 if (data.getIndices().size() != data.getN().size2()) {
1014 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1015 "data inconsistency");
1016 }
1017 int nb_row = data.getIndices().size();
1018 Nf.resize(nb_row);
1019 Nf.clear();
1020 int nb_gauss_pts = data.getN().size1();
1021 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1022 double w = getGaussPts()(3, gg) * getVolume();
1023 noalias(Nf) += w * data.getN(gg) * cTx.divergenceAtGaussPts[gg];
1024 }
1025 CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
1026 ;
1028 }
1029 };
1030
1031 /** \brief Calculate source therms, i.e. \f$\int_\mathcal{T} f v
1032 * \textrm{d}\mathcal{T}\f$
1033 */
1038
1039 OpL2Source(MixTransportElement &ctx, const std::string val_name, Vec f)
1040 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1041 val_name, UserDataOperator::OPROW),
1042 cTx(ctx), F(f) {}
1043
1048 if (data.getFieldData().size() == 0)
1050 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1051 int nb_row = data.getFieldData().size();
1052 Nf.resize(nb_row, false);
1053 Nf.clear();
1054 int nb_gauss_pts = data.getN().size1();
1055 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1056 double w = getGaussPts()(3, gg) * getVolume();
1057 const double x = getCoordsAtGaussPts()(gg, 0);
1058 const double y = getCoordsAtGaussPts()(gg, 1);
1059 const double z = getCoordsAtGaussPts()(gg, 2);
1060 double flux = 0;
1061 CHKERR cTx.getSource(fe_ent, x, y, z, flux);
1062 noalias(Nf) += w * data.getN(gg) * flux;
1063 }
1064 CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
1065 ;
1067 }
1068 };
1069
1070 /**
1071 * \brief calculate \f$ \int_\mathcal{S} {\boldsymbol\tau} \cdot \mathbf{n}u
1072 \textrm{d}\mathcal{S} \f$
1073
1074 * This terms comes from differentiation by parts. Note that in this Dirichlet
1075 * boundary conditions are natural.
1076
1077 */
1080
1083
1084 /**
1085 * \brief Constructor
1086 */
1087 OpRhsBcOnValues(MixTransportElement &ctx, const std::string fluxes_name,
1088 Vec f)
1090 fluxes_name, UserDataOperator::OPROW),
1091 cTx(ctx), F(f) {}
1092
1094
1095 /**
1096 * \brief Integrate boundary condition
1097 * @param side local index of entity
1098 * @param type type of entity
1099 * @param data data on entity
1100 * @return error code
1101 */
1105 if (data.getFieldData().size() == 0)
1107 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1108 nF.resize(data.getIndices().size());
1109 nF.clear();
1110 int nb_gauss_pts = data.getN().size1();
1111 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1112 const double x = getCoordsAtGaussPts()(gg, 0);
1113 const double y = getCoordsAtGaussPts()(gg, 1);
1114 const double z = getCoordsAtGaussPts()(gg, 2);
1115
1116 double value;
1117 CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
1118 ;
1119 double w = getGaussPts()(2, gg) * 0.5;
1120 if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1121 noalias(nF) +=
1122 w * prod(data.getVectorN<3>(gg), getNormalsAtGaussPts(gg)) *
1123 value;
1124 } else {
1125 noalias(nF) += w * prod(data.getVectorN<3>(gg), getNormal()) * value;
1126 }
1127 }
1128 Vec f = (F != PETSC_NULL) ? F : getFEMethod()->ts_F;
1129 CHKERR VecSetValues(f, data.getIndices().size(), &data.getIndices()[0],
1130 &nF[0], ADD_VALUES);
1131
1133 }
1134 };
1135
1136 /**
1137 * \brief Evaluate boundary conditions on fluxes.
1138 *
1139 * Note that Neumann boundary conditions here are essential. So it is opposite
1140 * what you find in displacement finite element method.
1141 *
1142
1143 * Here we have to solve for degrees of freedom on boundary such base
1144 functions
1145 * approximate flux.
1146
1147 *
1148 */
1152 OpEvaluateBcOnFluxes(MixTransportElement &ctx, const std::string flux_name)
1154 flux_name, UserDataOperator::OPROW),
1155 cTx(ctx) {}
1156
1160
1164 if (data.getFieldData().size() == 0)
1166 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1167 int nb_dofs = data.getFieldData().size();
1168 int nb_gauss_pts = data.getN().size1();
1169 if (3 * nb_dofs != static_cast<int>(data.getN().size2())) {
1170 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1171 "wrong number of dofs");
1172 }
1173 NN.resize(nb_dofs, nb_dofs);
1174 Nf.resize(nb_dofs);
1175 NN.clear();
1176 Nf.clear();
1177
1178 // Get normal vector. Note that when higher order geometry is set, then
1179 // face element could be curved, i.e. normal can be different at each
1180 // integration point.
1181 double *normal_ptr;
1182 if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1183 // HO geometry
1184 normal_ptr = &getNormalsAtGaussPts(0)[0];
1185 } else {
1186 // Linear geometry, i.e. constant normal on face
1187 normal_ptr = &getNormal()[0];
1188 }
1189 // set tensor from pointer
1190 FTensor::Tensor1<const double *, 3> t_normal(normal_ptr, &normal_ptr[1],
1191 &normal_ptr[2], 3);
1192 // get base functions
1193 auto t_n_hdiv_row = data.getFTensor1N<3>();
1194
1195 double nrm2 = 0;
1196
1197 // loop over integration points
1198 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1199
1200 // get integration point coordinates
1201 const double x = getCoordsAtGaussPts()(gg, 0);
1202 const double y = getCoordsAtGaussPts()(gg, 1);
1203 const double z = getCoordsAtGaussPts()(gg, 2);
1204
1205 // get flux on fece for given element handle and coordinates
1206 double flux;
1207 CHKERR cTx.getBcOnFluxes(fe_ent, x, y, z, flux);
1208 ;
1209 // get weight for integration rule
1210 double w = getGaussPts()(2, gg);
1211 if (gg == 0) {
1212 nrm2 = sqrt(t_normal(i) * t_normal(i));
1213 }
1214
1215 // set tensor of rank 0 to matrix NN elements
1216 // loop over base functions on rows and columns
1217 for (int ll = 0; ll != nb_dofs; ll++) {
1218 // get column on shape functions
1220 &data.getVectorN<3>(gg)(0, HVEC0),
1221 &data.getVectorN<3>(gg)(0, HVEC1),
1222 &data.getVectorN<3>(gg)(0, HVEC2), 3);
1223 for (int kk = 0; kk <= ll; kk++) {
1224 NN(ll, kk) += w * t_n_hdiv_row(i) * t_n_hdiv_col(i);
1225 ++t_n_hdiv_col;
1226 }
1227 // right hand side
1228 Nf[ll] += w * t_n_hdiv_row(i) * t_normal(i) * flux / nrm2;
1229 ++t_n_hdiv_row;
1230 }
1231
1232 // If HO geometry increment t_normal to next integration point
1233 if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1234 ++t_normal;
1235 nrm2 = sqrt(t_normal(i) * t_normal(i));
1236 }
1237 }
1238 // get global dofs indices on element
1239 cTx.bcIndices.insert(data.getIndices().begin(), data.getIndices().end());
1240 // factor matrix
1242 // solve local problem
1243 cholesky_solve(NN, Nf, ublas::lower());
1244
1245 // set solution to vector
1246 CHKERR VecSetOption(cTx.D0, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1247 CHKERR VecSetValues(cTx.D0, nb_dofs, &*data.getIndices().begin(),
1248 &*Nf.begin(), INSERT_VALUES);
1249 for (int dd = 0; dd != nb_dofs; ++dd)
1250 data.getFieldDofs()[dd]->getFieldData() = Nf[dd];
1251
1253 }
1254 };
1255
1256 /**
1257 * \brief Calculate values at integration points
1258 */
1261
1263
1265 const std::string val_name = "VALUES")
1266 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1267 val_name, UserDataOperator::OPROW),
1268 cTx(ctx) {}
1269
1271
1275 if (data.getFieldData().size() == 0)
1277
1278 int nb_gauss_pts = data.getN().size1();
1279 cTx.valuesAtGaussPts.resize(nb_gauss_pts);
1280 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1281 cTx.valuesAtGaussPts[gg] =
1282 inner_prod(trans(data.getN(gg)), data.getFieldData());
1283 }
1284
1286 }
1287 };
1288
1289 /**
1290 * \brief Calculate gradients of values at integration points
1291 */
1294
1296
1298 const std::string val_name = "VALUES")
1299 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1300 val_name, UserDataOperator::OPROW),
1301 cTx(ctx) {}
1303
1307 if (data.getFieldData().size() == 0)
1309 int nb_gauss_pts = data.getDiffN().size1();
1310 cTx.valuesGradientAtGaussPts.resize(3, nb_gauss_pts);
1311 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1312 ublas::matrix_column<MatrixDouble> values_grad_at_gauss_pts(
1313 cTx.valuesGradientAtGaussPts, gg);
1314 noalias(values_grad_at_gauss_pts) =
1315 prod(trans(data.getDiffN(gg)), data.getFieldData());
1316 }
1318 }
1319 };
1320
1321 /**
1322 * \brief calculate flux at integration point
1323 */
1326
1328
1330 const std::string field_name)
1331 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1333 cTx(ctx) {}
1334
1339
1340 if (data.getFieldData().size() == 0)
1342 int nb_gauss_pts = data.getDiffN().size1();
1343 int nb_dofs = data.getFieldData().size();
1344 cTx.fluxesAtGaussPts.resize(3, nb_gauss_pts);
1345 cTx.divergenceAtGaussPts.resize(nb_gauss_pts);
1346 if (type == MBTRI && side == 0) {
1347 cTx.divergenceAtGaussPts.clear();
1348 cTx.fluxesAtGaussPts.clear();
1349 }
1350 divVec.resize(nb_dofs);
1351
1353 auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
1354
1355 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1356 for(auto &v : divVec) {
1357 v = t_base_diff_hdiv(i, i);
1358 ++t_base_diff_hdiv;
1359 }
1360 cTx.divergenceAtGaussPts[gg] += inner_prod(divVec, data.getFieldData());
1361 ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1362 cTx.fluxesAtGaussPts, gg);
1363 flux_at_gauss_pt +=
1364 prod(trans(data.getVectorN<3>(gg)), data.getFieldData());
1365 }
1367 }
1368 };
1369
1370 map<double, EntityHandle> errorMap;
1374
1375 /** \brief calculate error evaluator
1376 */
1377 struct OpError
1379
1381
1382 OpError(MixTransportElement &ctx, const std::string field_name)
1383 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1385 cTx(ctx) {}
1386 virtual ~OpError() {}
1387
1390
1394 if (type != MBTET)
1396 invK.resize(3, 3, false);
1397 int nb_gauss_pts = data.getN().size1();
1398 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1399 double def_val = 0;
1400 Tag th_error_flux, th_error_div;
1401 CHKERR cTx.mField.get_moab().tag_get_handle(
1402 "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
1403 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1404 double *error_flux_ptr;
1405 CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1406 th_error_flux, &fe_ent, 1, (const void **)&error_flux_ptr);
1407
1408 CHKERR cTx.mField.get_moab().tag_get_handle(
1409 "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
1410 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1411 double *error_div_ptr;
1412 CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1413 th_error_div, &fe_ent, 1, (const void **)&error_div_ptr);
1414
1415 Tag th_error_jump;
1416 CHKERR cTx.mField.get_moab().tag_get_handle(
1417 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1418 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1419 double *error_jump_ptr;
1420 CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1421 th_error_jump, &fe_ent, 1, (const void **)&error_jump_ptr);
1422 *error_jump_ptr = 0;
1423
1424 /// characteristic size of the element
1425 const double h = pow(getVolume() * 12 / std::sqrt(2), (double)1 / 3);
1426
1427 for (int ff = 0; ff != 4; ff++) {
1428 EntityHandle face;
1429 CHKERR cTx.mField.get_moab().side_element(fe_ent, 2, ff, face);
1430 double *error_face_jump_ptr;
1431 CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1432 th_error_jump, &face, 1, (const void **)&error_face_jump_ptr);
1433 *error_face_jump_ptr = (1 / sqrt(h)) * sqrt(*error_face_jump_ptr);
1434 *error_face_jump_ptr = pow(*error_face_jump_ptr, 2);
1435 *error_jump_ptr += *error_face_jump_ptr;
1436 }
1437
1438 *error_flux_ptr = 0;
1439 *error_div_ptr = 0;
1440 deltaFlux.resize(3, false);
1441 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1442 double w = getGaussPts()(3, gg) * getVolume();
1443 const double x = getCoordsAtGaussPts()(gg, 0);
1444 const double y = getCoordsAtGaussPts()(gg, 1);
1445 const double z = getCoordsAtGaussPts()(gg, 2);
1446 double flux;
1447 CHKERR cTx.getSource(fe_ent, x, y, z, flux);
1448 ;
1449 *error_div_ptr += w * pow(cTx.divergenceAtGaussPts[gg] - flux, 2);
1450 CHKERR cTx.getResistivity(fe_ent, x, y, z, invK);
1451 ;
1452 ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1453 cTx.fluxesAtGaussPts, gg);
1454 ublas::matrix_column<MatrixDouble> values_grad_at_gauss_pts(
1455 cTx.valuesGradientAtGaussPts, gg);
1456 noalias(deltaFlux) =
1457 prod(invK, flux_at_gauss_pt) + values_grad_at_gauss_pts;
1458 *error_flux_ptr += w * inner_prod(deltaFlux, deltaFlux);
1459 }
1460 *error_div_ptr = h * sqrt(*error_div_ptr);
1461 *error_div_ptr = pow(*error_div_ptr, 2);
1462 cTx.errorMap[sqrt(*error_flux_ptr + *error_div_ptr + *error_jump_ptr)] =
1463 fe_ent;
1464 // Sum/Integrate all errors
1465 cTx.sumErrorFlux += *error_flux_ptr * getVolume();
1466 cTx.sumErrorDiv += *error_div_ptr * getVolume();
1467 // FIXME: Summation should be while skeleton is calculated
1468 cTx.sumErrorJump +=
1469 *error_jump_ptr * getVolume(); /// FIXME: this need to be fixed
1471 }
1472 };
1473
1474 /**
1475 * \brief calculate jump on entities
1476 */
1479
1480 /**
1481 * \brief volume element to get values from adjacent tets to face
1482 */
1483 VolumeElementForcesAndSourcesCoreOnSide volSideFe;
1484
1485 /** store values at integration point, key of the map is sense of face in
1486 * respect to adjacent tetrahedra
1487 */
1488 map<int, VectorDouble> valMap;
1489
1490 /**
1491 * \brief calculate values on adjacent tetrahedra to face
1492 */
1495 map<int, VectorDouble> &valMap;
1496 OpVolSide(map<int, VectorDouble> &val_map)
1497 : VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator(
1498 "VALUES", UserDataOperator::OPROW),
1499 valMap(val_map) {}
1503 if (data.getFieldData().size() == 0)
1505 int nb_gauss_pts = data.getN().size1();
1506 valMap[getFaceSense()].resize(nb_gauss_pts);
1507 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1508 valMap[getFaceSense()][gg] =
1509 inner_prod(trans(data.getN(gg)), data.getFieldData());
1510 }
1512 }
1513 };
1514
1516
1517 OpSkeleton(MixTransportElement &ctx, const std::string flux_name)
1519 flux_name, UserDataOperator::OPROW),
1520 volSideFe(ctx.mField), cTx(ctx) {
1521 volSideFe.getOpPtrVector().push_back(new OpSkeleton::OpVolSide(valMap));
1522 }
1523
1527
1528 if (type == MBTRI) {
1529 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1530
1531 double def_val = 0;
1532 Tag th_error_jump;
1533 CHKERR cTx.mField.get_moab().tag_get_handle(
1534 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1535 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1536 double *error_jump_ptr;
1537 CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1538 th_error_jump, &fe_ent, 1, (const void **)&error_jump_ptr);
1539 *error_jump_ptr = 0;
1540
1541 // check if this is essential boundary condition
1542 EntityHandle essential_bc_meshset =
1543 cTx.mField.get_finite_element_meshset("MIX_BCFLUX");
1544 if (cTx.mField.get_moab().contains_entities(essential_bc_meshset,
1545 &fe_ent, 1)) {
1546 // essential bc, np jump then, exit and go to next face
1548 }
1549
1550 // calculate values form adjacent tets
1551 valMap.clear();
1552 CHKERR loopSideVolumes("MIX", volSideFe);
1553 ;
1554
1555 int nb_gauss_pts = data.getN().size1();
1556
1557 // it is only one face, so it has to be bc natural boundary condition
1558 if (valMap.size() == 1) {
1559 if (static_cast<int>(valMap.begin()->second.size()) != nb_gauss_pts) {
1560 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1561 "wrong number of integration points");
1562 }
1563 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1564 const double x = getCoordsAtGaussPts()(gg, 0);
1565 const double y = getCoordsAtGaussPts()(gg, 1);
1566 const double z = getCoordsAtGaussPts()(gg, 2);
1567 double value;
1568 CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
1569 double w = getGaussPts()(2, gg);
1570 if (static_cast<int>(getNormalsAtGaussPts().size1()) ==
1571 nb_gauss_pts) {
1572 w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1573 } else {
1574 w *= getArea();
1575 }
1576 *error_jump_ptr += w * pow(value - valMap.begin()->second[gg], 2);
1577 }
1578 } else if (valMap.size() == 2) {
1579 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1580 double w = getGaussPts()(2, gg);
1581 if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1582 w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1583 } else {
1584 w *= getArea();
1585 }
1586 double delta = valMap.at(1)[gg] - valMap.at(-1)[gg];
1587 *error_jump_ptr += w * pow(delta, 2);
1588 }
1589 } else {
1590 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1591 "data inconsistency, wrong number of neighbors "
1592 "valMap.size() = %d",
1593 valMap.size());
1594 }
1595 }
1596
1598 }
1599 };
1600};
1601
1602} // namespace MixTransport
1603
1604#endif //_MIX_TRANPORT_ELEMENT_HPP_
1605
1606/***************************************************************************/ /**
1607 * \defgroup mofem_mix_transport_elem Mix transport element
1608 * \ingroup user_modules
1609 ******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
static const double eps
EntitiesFieldData::EntData EntData
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MAT_THERMALSET
block name is "MAT_THERMAL"
Definition: definitions.h:161
@ BLOCKSET
Definition: definitions.h:148
@ HVEC0
Definition: definitions.h:186
@ HVEC1
Definition: definitions.h:186
@ HVEC2
Definition: definitions.h:186
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual EntityHandle get_finite_element_meshset(const std::string &name) const =0
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
auto f
Definition: HenckyOps.hpp:5
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
double w(const double g, const double t)
const double D
diffusivity
double h
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
static constexpr double delta
data for evaluation of het conductivity and heat capacity elements
Range tEts
constatins elements in block set
OpDivTauU_HdivL2(MixTransportElement &ctx, const std::string flux_name_row, string val_name_col, Vec f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpError(MixTransportElement &ctx, const std::string field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpEvaluateBcOnFluxes(MixTransportElement &ctx, const std::string flux_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpFluxDivergenceAtGaussPts(MixTransportElement &ctx, const std::string field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpL2Source(MixTransportElement &ctx, const std::string val_name, Vec f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpPostProc(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts)
OpRhsBcOnValues(MixTransportElement &ctx, const std::string fluxes_name, Vec f)
Constructor.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
calculate values on adjacent tetrahedra to face
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpSkeleton(MixTransportElement &ctx, const std::string flux_name)
VolumeElementForcesAndSourcesCoreOnSide volSideFe
volume element to get values from adjacent tets to face
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
OpTauDotSigma_HdivHdiv(MixTransportElement &ctx, const std::string flux_name, Mat aij, Vec f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Assemble matrix.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
OpVDivSigma_L2Hdiv(MixTransportElement &ctx, const std::string val_name_row, string flux_name_col, Mat aij, Vec f)
Constructor.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpValuesAtGaussPts(MixTransportElement &ctx, const std::string val_name="VALUES")
Calculate gradients of values at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpValuesGradientAtGaussPts(MixTransportElement &ctx, const std::string val_name="VALUES")
MoFEMErrorCode postProc(const string out_file)
Post process results.
MoFEMErrorCode destroyMatrices()
destroy matrices
MoFEMErrorCode evaluateError()
Calculate error on elements.
MoFEMErrorCode solveLinearProblem()
solve problem
MoFEMErrorCode getDirichletBCIndices(IS *is)
get dof indices where essential boundary conditions are applied
MoFEMErrorCode createMatrices()
create matrices
MatrixDouble valuesGradientAtGaussPts
gradients at integration points on element
virtual MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
MoFEMErrorCode calculateResidual()
calculate residual
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
VectorDouble valuesAtGaussPts
values at integration points on element
map< double, EntityHandle > errorMap
MyTriFE feTri
Instance of surface element.
virtual MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
evaluate natural (Dirichlet) boundary conditions
MyVolumeFE feVol
Instance of volume element.
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
Add fields to database.
virtual MoFEMErrorCode getResistivity(const EntityHandle ent, const double x, const double y, const double z, MatrixDouble3by3 &inv_k)
natural (Dirichlet) boundary conditions (set values)
MixTransportElement(MoFEM::Interface &m_field)
construction of this data structure
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
VectorDouble divergenceAtGaussPts
divergence at integration points on element
virtual MoFEMErrorCode getSource(const EntityHandle ent, const double x, const double y, const double z, double &flux)
set source term
MatrixDouble fluxesAtGaussPts
fluxes at integration points on element
MoFEMErrorCode buildProblem(BitRefLevel &ref_level)
Build problem.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Deprecated interface functions.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
@ OPCOL
operator doWork function is executed on FE columns
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Post post-proc data at points from hash maps.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCore * getVolumeFE() const
return pointer to Generic Volume Finite Element object
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)