v0.14.0
Searching...
No Matches
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#ifndef _MIX_TRANPORT_ELEMENT_HPP_
9#define _MIX_TRANPORT_ELEMENT_HPP_
10
11namespace MixTransport {
12
13/** \brief Mix transport problem
14 \ingroup mofem_mix_transport_elem
15
16 Note to solve this system you need to use direct solver or proper
18
19 It is based on \cite arnold2006differential \cite arnold2012mixed \cite
20 reviartthomas1996
22
23 General problem have form,
24 \f[
25 \mathbf{A} \boldsymbol\sigma + \textrm{grad}[u] = \mathbf{0} \; \textrm{on} \;
26 \Omega \\ \textrm{div}[\boldsymbol\sigma] = f \; \textrm{on} \; \Omega \f]
27
28*/
30
32
33 /**
34 * \brief definition of volume element
35
36 * It is used to calculate volume integrals. On volume element we set-up
37 * operators to calculate components of matrix and vector.
38
39 */
43 int getRule(int order) { return 2 * order + 1; };
44 };
45
46 MyVolumeFE feVol; ///< Instance of volume element
47
48 /** \brief definition of surface element
49
50 * It is used to calculate surface integrals. On volume element are operators
51 * evaluating natural boundary conditions.
52
53 */
57 int getRule(int order) { return 2 * order + 1; };
58 };
59
60 MyTriFE feTri; ///< Instance of surface element
61
62 /**
63 * \brief construction of this data structure
64 */
66 : mField(m_field), feVol(m_field), feTri(m_field){};
67
68 /**
69 * \brief destructor
70 */
72
73 VectorDouble valuesAtGaussPts; ///< values at integration points on element
74 MatrixDouble
76 VectorDouble
77 divergenceAtGaussPts; ///< divergence at integration points on element
78 MatrixDouble fluxesAtGaussPts; ///< fluxes at integration points on element
79
80 set<int> bcIndices;
81
82 /**
83 * \brief get dof indices where essential boundary conditions are applied
84 * @param is indices
85 * @return error code
86 */
87 MoFEMErrorCode getDirichletBCIndices(IS *is) {
89 std::vector<int> ids;
90 ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
91 IS is_local;
92 CHKERR ISCreateGeneral(mField.get_comm(), ids.size(),
93 ids.empty() ? PETSC_NULL : &ids[0],
94 PETSC_COPY_VALUES, &is_local);
95 CHKERR ISAllGather(is_local, is);
96 CHKERR ISDestroy(&is_local);
98 }
99
100 /**
101 * \brief set source term
102 * @param ent handle to entity on which function is evaluated
103 * @param x coord
104 * @param y coord
105 * @param z coord
106 * @param flux reference to source term set by function
107 * @return error code
108 */
109 virtual MoFEMErrorCode getSource(const EntityHandle ent, const double x,
110 const double y, const double z,
111 double &flux) {
113 flux = 0;
115 }
116
117 /**
118 * \brief natural (Dirichlet) boundary conditions (set values)
119 * @param ent handle to finite element entity
120 * @param x coord
121 * @param y coord
122 * @param z coord
123 * @param value reference to value set by function
124 * @return error code
125 */
126 virtual MoFEMErrorCode getResistivity(const EntityHandle ent, const double x,
127 const double y, const double z,
128 MatrixDouble3by3 &inv_k) {
130 inv_k.clear();
131 for (int dd = 0; dd < 3; dd++) {
132 inv_k(dd, dd) = 1;
133 }
135 }
136
137 /**
138 * \brief evaluate natural (Dirichlet) boundary conditions
139 * @param ent entity on which bc is evaluated
140 * @param x coordinate
141 * @param y coordinate
142 * @param z coordinate
143 * @param value vale
144 * @return error code
145 */
146 virtual MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg,
147 const double x, const double y,
148 const double z, double &value) {
150 value = 0;
152 }
153
154 /**
155 * \brief essential (Neumann) boundary condition (set fluxes)
156 * @param ent handle to finite element entity
157 * @param x coord
158 * @param y coord
159 * @param z coord
160 * @param flux reference to flux which is set by function
161 * @return [description]
162 */
163 virtual MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x,
164 const double y, const double z,
165 double &flux) {
167 flux = 0;
169 }
170
171 /** \brief data for evaluation of het conductivity and heat capacity elements
172 * \ingroup mofem_thermal_elem
173 */
174 struct BlockData {
176 double cApacity;
177 Range tEts; ///< constatins elements in block set
178 };
179 std::map<int, BlockData>
180 setOfBlocks; ///< maps block set id with appropriate BlockData
181
182 /**
183 * \brief Add fields to database
184 * @param values name of the fields
185 * @param fluxes name of filed for fluxes
186 * @param order order of approximation
187 * @return error code
188 */
189 MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes,
190 const int order) {
192 // Fields
195
196 // meshset consisting all entities in mesh
197 EntityHandle root_set = mField.get_moab().get_root_set();
198 // add entities to field
201 CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
202 CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
203 CHKERR mField.set_field_order(root_set, MBTET, values, order);
205 }
206
207 /// \brief add finite elements
209 const std::string &values_name) {
211
212 // Set up volume element operators. Operators are used to calculate
213 // components of stiffness matrix & right hand side, in essence are used to
214 // do volume integrals over tetrahedral in this case.
215
216 // Define element "MIX". Note that this element will work with fluxes_name
217 // and values_name. This reflect bilinear form for the problem
225
226 // Define finite element to integrate over skeleton, we need that to
227 // evaluate error
230 fluxes_name);
232 fluxes_name);
234 fluxes_name);
235
236 // Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those
237 // BLOCKSETS and add them to "MIX" finite element. In addition get data form
238 // that meshset and set conductivity which is used to calculate fluxes from
239 // gradients of concentration or gradient of temperature, depending how you
240 // interpret variables.
243
244 Mat_Thermal temp_data;
245 CHKERR it->getAttributeDataStructure(temp_data);
246 setOfBlocks[it->getMeshsetId()].cOnductivity =
247 temp_data.data.Conductivity;
248 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
249 CHKERR mField.get_moab().get_entities_by_type(
250 it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
252 setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
253
254 Range skeleton;
256 setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
257 moab::Interface::UNION);
259 "MIX_SKELETON");
260 }
261
262 // Define element to integrate natural boundary conditions, i.e. set values.
265 fluxes_name);
267 fluxes_name);
269 fluxes_name);
271 values_name);
272
273 // Define element to apply essential boundary conditions.
276 fluxes_name);
278 fluxes_name);
280 fluxes_name);
282 values_name);
284 }
285
286 /**
287 * \brief Build problem
288 * @param ref_level mesh refinement on which mesh problem you like to built.
289 * @return error code
290 */
291 MoFEMErrorCode buildProblem(BitRefLevel &ref_level) {
293 // build field
295 // get tetrahedrons which has been build previously and now in so called
296 // garbage bit level
297 Range done_tets;
298 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
299 BitRefLevel().set(0), BitRefLevel().set(), MBTET, done_tets);
300 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
301 BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTET,
302 done_tets);
303 // get tetrahedrons which belong to problem bit level
304 Range ref_tets;
305 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
306 ref_level, BitRefLevel().set(), MBTET, ref_tets);
307 ref_tets = subtract(ref_tets, done_tets);
308 CHKERR mField.build_finite_elements("MIX", &ref_tets, 2);
309 // get triangles which has been build previously and now in so called
310 // garbage bit level
311 Range done_faces;
312 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
313 BitRefLevel().set(0), BitRefLevel().set(), MBTRI, done_faces);
314 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
315 BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTRI,
316 done_faces);
317 // get triangles which belong to problem bit level
318 Range ref_faces;
319 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
320 ref_level, BitRefLevel().set(), MBTRI, ref_faces);
321 ref_faces = subtract(ref_faces, done_faces);
322 // build finite elements structures
323 CHKERR mField.build_finite_elements("MIX_BCFLUX", &ref_faces, 2);
324 CHKERR mField.build_finite_elements("MIX_BCVALUE", &ref_faces, 2);
325 CHKERR mField.build_finite_elements("MIX_SKELETON", &ref_faces, 2);
326 // Build adjacencies of degrees of freedom and elements
328 // Define problem
330 // set refinement level for problem
332 // Add element to problem
335 // Boundary conditions
338 // build problem
339
340 ProblemsManager *prb_mng_ptr;
341 CHKERR mField.getInterface(prb_mng_ptr);
342 CHKERR prb_mng_ptr->buildProblem("MIX", true);
343 // mesh partitioning
344 // partition
345 CHKERR prb_mng_ptr->partitionProblem("MIX");
346 CHKERR prb_mng_ptr->partitionFiniteElements("MIX");
347 // what are ghost nodes, see Petsc Manual
348 CHKERR prb_mng_ptr->partitionGhostDofs("MIX");
350 }
351
354 moab::Interface &postProcMesh;
355 std::vector<EntityHandle> &mapGaussPts;
356 OpPostProc(moab::Interface &post_proc_mesh,
357 std::vector<EntityHandle> &map_gauss_pts)
358 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
359 "VALUES", UserDataOperator::OPCOL),
360 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
361 MoFEMErrorCode doWork(int side, EntityType type,
362 EntitiesFieldData::EntData &data) {
364 if (type != MBTET)
366 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
367 Tag th_error_flux;
368 CHKERR getVolumeFE()->mField.get_moab().tag_get_handle("ERROR_FLUX",
369 th_error_flux);
370 double *error_flux_ptr;
371 CHKERR getVolumeFE()->mField.get_moab().tag_get_by_ptr(
372 th_error_flux, &fe_ent, 1, (const void **)&error_flux_ptr);
373
374 Tag th_error_div;
375 CHKERR getVolumeFE()->mField.get_moab().tag_get_handle("ERROR_DIV",
376 th_error_div);
377 double *error_div_ptr;
378 CHKERR getVolumeFE()->mField.get_moab().tag_get_by_ptr(
379 th_error_div, &fe_ent, 1, (const void **)&error_div_ptr);
380
381 Tag th_error_jump;
382 CHKERR getVolumeFE()->mField.get_moab().tag_get_handle("ERROR_JUMP",
383 th_error_jump);
384 double *error_jump_ptr;
385 CHKERR getVolumeFE()->mField.get_moab().tag_get_by_ptr(
386 th_error_jump, &fe_ent, 1, (const void **)&error_jump_ptr);
387
388 {
389 double def_val = 0;
390 Tag th_error_flux;
391 CHKERR postProcMesh.tag_get_handle(
392 "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
393 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
394 for (vector<EntityHandle>::iterator vit = mapGaussPts.begin();
395 vit != mapGaussPts.end(); vit++) {
396 CHKERR postProcMesh.tag_set_data(th_error_flux, &*vit, 1,
397 error_flux_ptr);
398 }
399
400 Tag th_error_div;
401 CHKERR postProcMesh.tag_get_handle(
402 "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
403 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
404 for (vector<EntityHandle>::iterator vit = mapGaussPts.begin();
405 vit != mapGaussPts.end(); vit++) {
406 CHKERR postProcMesh.tag_set_data(th_error_div, &*vit, 1,
407 error_div_ptr);
408 }
409
410 Tag th_error_jump;
411 CHKERR postProcMesh.tag_get_handle(
412 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
413 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
414 for (vector<EntityHandle>::iterator vit = mapGaussPts.begin();
415 vit != mapGaussPts.end(); vit++) {
416 CHKERR postProcMesh.tag_set_data(th_error_jump, &*vit, 1,
417 error_jump_ptr);
418 }
419 }
421 }
422 };
423
424 /**
425 * \brief Post process results
426 * @return error code
427 */
428 MoFEMErrorCode postProc(const string out_file) {
430
431 PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
432 mField);
433
435
436 auto values_ptr = boost::make_shared<VectorDouble>();
438 auto flux_ptr = boost::make_shared<MatrixDouble>();
439
440 post_proc.getOpPtrVector().push_back(
441 new OpCalculateScalarFieldValues("VALUES", values_ptr));
442 post_proc.getOpPtrVector().push_back(
444 post_proc.getOpPtrVector().push_back(
445 new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
446
447 using OpPPMap = OpPostProcMapInMoab<3, 3>;
448
449 post_proc.getOpPtrVector().push_back(
450
451 new OpPPMap(post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
452
453 {{"VALUES", values_ptr}},
454
456
457 {}, {}
458
459 ));
460
461 post_proc.getOpPtrVector().push_back(new OpPostProc(
462 post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
463
464 CHKERR mField.loop_finite_elements("MIX", "MIX", post_proc);
465
466 CHKERR post_proc.writeFile(out_file.c_str());
468 }
469
470 Vec D, D0, F;
471 Mat Aij;
472
473 /// \brief create matrices
474 MoFEMErrorCode createMatrices() {
476 CHKERR mField.getInterface<MatrixManager>()
477 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("MIX", &Aij);
478 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
479 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
480 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
482 }
483
484 /**
485 * \brief solve problem
486 * @return error code
487 */
488 MoFEMErrorCode solveLinearProblem() {
490 CHKERR MatZeroEntries(Aij);
491 CHKERR VecZeroEntries(F);
492 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
493 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
494 CHKERR VecZeroEntries(D0);
495 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
496 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
497 CHKERR VecZeroEntries(D);
498 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
499 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
500
501 CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
502 "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
503
504 // Calculate essential boundary conditions
505
506 // clear essential bc indices, it could have dofs from other mesh refinement
507 bcIndices.clear();
508 // clear operator, just in case if some other operators are left on this
509 // element
510 feTri.getOpPtrVector().clear();
512 // set operator to calculate essential boundary conditions
513 feTri.getOpPtrVector().push_back(new OpEvaluateBcOnFluxes(*this, "FLUXES"));
514 CHKERR mField.loop_finite_elements("MIX", "MIX_BCFLUX", feTri);
515 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
516 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
517 CHKERR VecAssemblyBegin(D0);
518 CHKERR VecAssemblyEnd(D0);
519
520 // set operators to calculate matrix and right hand side vectors
521 feVol.getOpPtrVector().clear();
523 feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
524 feVol.getOpPtrVector().push_back(
525 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
526 feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
527 feVol.getOpPtrVector().push_back(
528 new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
529 feVol.getOpPtrVector().push_back(
530 new OpTauDotSigma_HdivHdiv(*this, "FLUXES", Aij, F));
531 feVol.getOpPtrVector().push_back(
532 new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", Aij, F));
533 CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
534
535 // calculate right hand side for natural boundary conditions
536 feTri.getOpPtrVector().clear();
538 feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
539 CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
540
541 // assemble matrices
542 CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
543 CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
546 CHKERR VecAssemblyBegin(F);
547 CHKERR VecAssemblyEnd(F);
548
549 {
550 double nrm2_F;
551 CHKERR VecNorm(F, NORM_2, &nrm2_F);
552 PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
553 }
554
556
557 // for ksp solver vector is moved into rhs side
558 // for snes it is left ond the left
559 CHKERR VecScale(F, -1);
560
561 IS essential_bc_ids;
562 CHKERR getDirichletBCIndices(&essential_bc_ids);
563 CHKERR MatZeroRowsColumnsIS(Aij, essential_bc_ids, 1, D0, F);
564 CHKERR ISDestroy(&essential_bc_ids);
565
566 // {
567 // double norm;
568 // CHKERR MatNorm(Aij,NORM_FROBENIUS,&norm); ;
569 // PetscPrintf(PETSC_COMM_WORLD,"mat norm = %6.4e\n",norm);
570 // }
571
572 {
573 double nrm2_F;
574 CHKERR VecNorm(F, NORM_2, &nrm2_F);
575 PetscPrintf(PETSC_COMM_WORLD, "With BC nrm2_F = %6.4e\n", nrm2_F);
576 }
577
578 // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
579 // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
580 // std::string wait;
581 // std::cin >> wait;
582
583 // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
584 // std::cin >> wait;
585
586 // Solve
587 KSP solver;
588 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
589 CHKERR KSPSetOperators(solver, Aij, Aij);
590 CHKERR KSPSetFromOptions(solver);
591 CHKERR KSPSetUp(solver);
592 CHKERR KSPSolve(solver, F, D);
593 CHKERR KSPDestroy(&solver);
594
595 {
596 double nrm2_D;
597 CHKERR VecNorm(D, NORM_2, &nrm2_D);
598 ;
599 PetscPrintf(PETSC_COMM_WORLD, "nrm2_D = %6.4e\n", nrm2_D);
600 }
601 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
602 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
603
604 // copy data form vector on mesh
605 CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
606 "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
607
609 }
610
611 /// \brief calculate residual
612 MoFEMErrorCode calculateResidual() {
614 CHKERR VecZeroEntries(F);
615 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
616 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
617 CHKERR VecAssemblyBegin(F);
618 CHKERR VecAssemblyEnd(F);
619 // calculate residuals
620 feVol.getOpPtrVector().clear();
622 feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
623 feVol.getOpPtrVector().push_back(
624 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
625 feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
626 feVol.getOpPtrVector().push_back(
627 new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
628 feVol.getOpPtrVector().push_back(
629 new OpTauDotSigma_HdivHdiv(*this, "FLUXES", PETSC_NULL, F));
630 feVol.getOpPtrVector().push_back(
631 new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", PETSC_NULL, F));
632 CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
633 feTri.getOpPtrVector().clear();
634 feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
635 CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
638 CHKERR VecAssemblyBegin(F);
639 CHKERR VecAssemblyEnd(F);
640 // CHKERR VecAXPY(F,-1.,D0);
641 // CHKERR MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULL,F);
642 {
643 std::vector<int> ids;
644 ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
645 std::vector<double> vals(ids.size(), 0);
646 CHKERR VecSetValues(F, ids.size(), &*ids.begin(), &*vals.begin(),
647 INSERT_VALUES);
648 CHKERR VecAssemblyBegin(F);
649 CHKERR VecAssemblyEnd(F);
650 }
651 {
652 double nrm2_F;
653 CHKERR VecNorm(F, NORM_2, &nrm2_F);
654 PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
655 const double eps = 1e-8;
656 if (nrm2_F > eps) {
657 // SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"problem with
658 // residual");
659 }
660 }
662 }
663
664 /** \brief Calculate error on elements
665
666 \todo this functions runs serial, in future should be parallel and work on
667 distributed meshes.
668
669 */
670 MoFEMErrorCode evaluateError() {
672 errorMap.clear();
673 sumErrorFlux = 0;
674 sumErrorDiv = 0;
675 sumErrorJump = 0;
676 feTri.getOpPtrVector().clear();
677 feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
678 CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
679 mField.get_comm_size());
680 feVol.getOpPtrVector().clear();
682 feVol.getOpPtrVector().push_back(
683 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
684 feVol.getOpPtrVector().push_back(
686 feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
687 CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
688 mField.get_comm_size());
689 const Problem *problem_ptr;
690 CHKERR mField.get_problem("MIX", &problem_ptr);
691 PetscPrintf(mField.get_comm(),
692 "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
693 "jump^2 = %6.4e error tot^2 = %6.4e\n",
694 problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
695 sumErrorJump, sumErrorFlux + sumErrorDiv + sumErrorJump);
697 }
698
699 /// \brief destroy matrices
700 MoFEMErrorCode destroyMatrices() {
702 CHKERR MatDestroy(&Aij);
703 ;
704 CHKERR VecDestroy(&D);
705 ;
706 CHKERR VecDestroy(&D0);
707 ;
708 CHKERR VecDestroy(&F);
709 ;
711 }
712
713 /**
714 \brief Assemble \f$\int_\mathcal{T} \mathbf{A} \boldsymbol\sigma \cdot 715 \boldsymbol\tau \textrm{d}\mathcal{T}\f$
716
717 \ingroup mofem_mix_transport_elem
718 */
721
723 Mat Aij;
724 Vec F;
725
727 const std::string flux_name, Mat aij, Vec f)
728 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
729 flux_name, flux_name,
730 UserDataOperator::OPROWCOL | UserDataOperator::OPCOL),
731 cTx(ctx), Aij(aij), F(f) {
732 sYmm = true;
733 }
734
735 MatrixDouble NN, transNN;
736 MatrixDouble3by3 invK;
737 VectorDouble Nf;
738
739 /**
740 * \brief Assemble matrix
741 * @param row_side local index of row entity on element
742 * @param col_side local index of col entity on element
743 * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
744 * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
745 * @param row_data data for row
746 * @param col_data data for col
747 * @return error code
748 */
749 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
750 EntityType col_type,
751 EntitiesFieldData::EntData &row_data,
752 EntitiesFieldData::EntData &col_data) {
754 if (Aij == PETSC_NULL)
756 if (row_data.getIndices().size() == 0)
758 if (col_data.getIndices().size() == 0)
760 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
761 int nb_row = row_data.getIndices().size();
762 int nb_col = col_data.getIndices().size();
763 NN.resize(nb_row, nb_col, false);
764 NN.clear();
767 invK.resize(3, 3, false);
770 &invK(0, 0), &invK(0, 1), &invK(0, 2), &invK(1, 0), &invK(1, 1),
771 &invK(1, 2), &invK(2, 0), &invK(2, 1), &invK(2, 2));
772 // Get base functions
773 auto t_n_hdiv_row = row_data.getFTensor1N<3>();
775 int nb_gauss_pts = row_data.getN().size1();
776 for (int gg = 0; gg != nb_gauss_pts; gg++) {
777 // get integration weight and multiply by element volume
778 double w = getGaussPts()(3, gg) * getVolume();
779 const double x = getCoordsAtGaussPts()(gg, 0);
780 const double y = getCoordsAtGaussPts()(gg, 1);
781 const double z = getCoordsAtGaussPts()(gg, 2);
782 // calculate receptivity (invers of conductivity)
783 CHKERR cTx.getResistivity(fe_ent, x, y, z, invK);
784 for (int kk = 0; kk != nb_row; kk++) {
786 &col_data.getVectorN<3>(gg)(0, HVEC0),
787 &col_data.getVectorN<3>(gg)(0, HVEC1),
788 &col_data.getVectorN<3>(gg)(0, HVEC2), 3);
789 t_row(j) = w * t_n_hdiv_row(i) * t_inv_k(i, j);
790 for (int ll = 0; ll != nb_col; ll++) {
791 NN(kk, ll) += t_row(j) * t_n_hdiv_col(j);
792 ++t_n_hdiv_col;
793 }
794 ++t_n_hdiv_row;
795 }
796 }
797 Mat a = (Aij != PETSC_NULL) ? Aij : getFEMethod()->ts_B;
798 CHKERR MatSetValues(a, nb_row, &row_data.getIndices()[0], nb_col,
800 // matrix is symmetric, assemble other part
801 if (row_side != col_side || row_type != col_type) {
802 transNN.resize(nb_col, nb_row);
803 noalias(transNN) = trans(NN);
804 CHKERR MatSetValues(a, nb_col, &col_data.getIndices()[0], nb_row,
805 &row_data.getIndices()[0], &transNN(0, 0),
807 }
808
810 }
811
812 /**
813 * \brief Assemble matrix
814 * @param side local index of row entity on element
815 * @param type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
816 * @param data data for row
817 * @return error code
818 */
819 MoFEMErrorCode doWork(int side, EntityType type,
820 EntitiesFieldData::EntData &data) {
822 if (F == PETSC_NULL)
824 int nb_row = data.getIndices().size();
825 if (nb_row == 0)
827
828 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
829 Nf.resize(nb_row);
830 Nf.clear();
833 invK.resize(3, 3, false);
834 Nf.resize(nb_row);
835 Nf.clear();
838 &invK(0, 0), &invK(0, 1), &invK(0, 2), &invK(1, 0), &invK(1, 1),
839 &invK(1, 2), &invK(2, 0), &invK(2, 1), &invK(2, 2));
840 // get base functions
841 auto t_n_hdiv = data.getFTensor1N<3>();
842 auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
843 int nb_gauss_pts = data.getN().size1();
844 for (int gg = 0; gg < nb_gauss_pts; gg++) {
845 double w = getGaussPts()(3, gg) * getVolume();
846 const double x = getCoordsAtGaussPts()(gg, 0);
847 const double y = getCoordsAtGaussPts()(gg, 1);
848 const double z = getCoordsAtGaussPts()(gg, 2);
849 CHKERR cTx.getResistivity(fe_ent, x, y, z, invK);
850 for (int ll = 0; ll != nb_row; ll++) {
851 Nf[ll] += w * t_n_hdiv(i) * t_inv_k(i, j) * t_flux(j);
852 ++t_n_hdiv;
853 }
854 ++t_flux;
855 }
856 CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
857
859 }
860 };
861
862 /** \brief Assemble \f$\int_\mathcal{T} u \textrm{div}[\boldsymbol\tau] 863 * \textrm{d}\mathcal{T} \f$
864 */
867
869 Vec F;
870
871 OpDivTauU_HdivL2(MixTransportElement &ctx, const std::string flux_name_row,
872 string val_name_col, Vec f)
873 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
874 flux_name_row, val_name_col, UserDataOperator::OPROW),
875 cTx(ctx), F(f) {
876 // this operator is not symmetric setting this variable makes element
877 // operator to loop over element entities (sub-simplices) without
878 // assumption that off-diagonal matrices are symmetric.
879 sYmm = false;
880 }
881
882 VectorDouble divVec, Nf;
883
884 MoFEMErrorCode doWork(int side, EntityType type,
885 EntitiesFieldData::EntData &data) {
887 if (data.getFieldData().size() == 0)
889 int nb_row = data.getIndices().size();
890 Nf.resize(nb_row);
891 Nf.clear();
892 divVec.resize(data.getN().size2() / 3, 0);
893 if (divVec.size() != data.getIndices().size()) {
894 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
895 "data inconsistency");
896 }
897 int nb_gauss_pts = data.getN().size1();
898
900 auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
901
902 int gg = 0;
903 for (; gg < nb_gauss_pts; gg++) {
904 double w = getGaussPts()(3, gg) * getVolume();
905 for (auto &v : divVec) {
906 v = t_base_diff_hdiv(i, i);
907 ++t_base_diff_hdiv;
908 }
909 noalias(Nf) -= w * divVec * cTx.valuesAtGaussPts[gg];
910 }
911 CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
912 ;
913
915 }
916 };
917
918 /** \brief \f$\int_\mathcal{T} \textrm{div}[\boldsymbol\sigma] v 919 * \textrm{d}\mathcal{T} \f$
920 */
923
925 Mat Aij;
926 Vec F;
927
928 /**
929 * \brief Constructor
930 */
931 OpVDivSigma_L2Hdiv(MixTransportElement &ctx, const std::string val_name_row,
932 string flux_name_col, Mat aij, Vec f)
933 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
934 val_name_row, flux_name_col,
935 UserDataOperator::OPROW | UserDataOperator::OPROWCOL),
936 cTx(ctx), Aij(aij), F(f) {
937
938 // this operator is not symmetric setting this variable makes element
939 // operator to loop over element entities without
940 // assumption that off-diagonal matrices are symmetric.
941 sYmm = false;
942 }
943
944 MatrixDouble NN, transNN;
945 VectorDouble divVec, Nf;
946
947 /**
948 * \brief Do calculations
949 * @param row_side local index of entity on row
950 * @param col_side local index of entity on column
951 * @param row_type type of row entity
952 * @param col_type type of col entity
953 * @param row_data row data structure carrying information about base
954 * functions, DOFs indices, etc.
955 * @param col_data column data structure carrying information about base
956 * functions, DOFs indices, etc.
957 * @return error code
958 */
959 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
960 EntityType col_type,
961 EntitiesFieldData::EntData &row_data,
962 EntitiesFieldData::EntData &col_data) {
964 if (Aij == PETSC_NULL)
966 if (row_data.getFieldData().size() == 0)
968 if (col_data.getFieldData().size() == 0)
970 int nb_row = row_data.getFieldData().size();
971 int nb_col = col_data.getFieldData().size();
972 NN.resize(nb_row, nb_col);
973 NN.clear();
974 divVec.resize(nb_col, false);
975
977 auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
978
979 int nb_gauss_pts = row_data.getN().size1();
980 for (int gg = 0; gg < nb_gauss_pts; gg++) {
981 double w = getGaussPts()(3, gg) * getVolume();
982 for (auto &v : divVec) {
983 v = t_base_diff_hdiv(i, i);
984 ++t_base_diff_hdiv;
985 }
986 noalias(NN) += w * outer_prod(row_data.getN(gg), divVec);
987 }
988 CHKERR MatSetValues(Aij, nb_row, &row_data.getIndices()[0], nb_col,
990 transNN.resize(nb_col, nb_row);
991 ublas::noalias(transNN) = -trans(NN);
992 CHKERR MatSetValues(Aij, nb_col, &col_data.getIndices()[0], nb_row,
993 &row_data.getIndices()[0], &transNN(0, 0),
996 }
997
998 MoFEMErrorCode doWork(int side, EntityType type,
999 EntitiesFieldData::EntData &data) {
1001 if (data.getIndices().size() == 0)
1003 if (data.getIndices().size() != data.getN().size2()) {
1004 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1005 "data inconsistency");
1006 }
1007 int nb_row = data.getIndices().size();
1008 Nf.resize(nb_row);
1009 Nf.clear();
1010 int nb_gauss_pts = data.getN().size1();
1011 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1012 double w = getGaussPts()(3, gg) * getVolume();
1013 noalias(Nf) += w * data.getN(gg) * cTx.divergenceAtGaussPts[gg];
1014 }
1015 CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
1016 ;
1018 }
1019 };
1020
1021 /** \brief Calculate source therms, i.e. \f$\int_\mathcal{T} f v 1022 * \textrm{d}\mathcal{T}\f$
1023 */
1027 Vec F;
1028
1029 OpL2Source(MixTransportElement &ctx, const std::string val_name, Vec f)
1030 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1031 val_name, UserDataOperator::OPROW),
1032 cTx(ctx), F(f) {}
1033
1034 VectorDouble Nf;
1035 MoFEMErrorCode doWork(int side, EntityType type,
1036 EntitiesFieldData::EntData &data) {
1038 if (data.getFieldData().size() == 0)
1040 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1041 int nb_row = data.getFieldData().size();
1042 Nf.resize(nb_row, false);
1043 Nf.clear();
1044 int nb_gauss_pts = data.getN().size1();
1045 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1046 double w = getGaussPts()(3, gg) * getVolume();
1047 const double x = getCoordsAtGaussPts()(gg, 0);
1048 const double y = getCoordsAtGaussPts()(gg, 1);
1049 const double z = getCoordsAtGaussPts()(gg, 2);
1050 double flux = 0;
1051 CHKERR cTx.getSource(fe_ent, x, y, z, flux);
1052 noalias(Nf) += w * data.getN(gg) * flux;
1053 }
1054 CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
1055 ;
1057 }
1058 };
1059
1060 /**
1061 * \brief calculate \f$\int_\mathcal{S} {\boldsymbol\tau} \cdot \mathbf{n}u 1062 \textrm{d}\mathcal{S} \f$
1063
1064 * This terms comes from differentiation by parts. Note that in this Dirichlet
1065 * boundary conditions are natural.
1066
1067 */
1070
1072 Vec F;
1073
1074 /**
1075 * \brief Constructor
1076 */
1077 OpRhsBcOnValues(MixTransportElement &ctx, const std::string fluxes_name,
1078 Vec f)
1080 fluxes_name, UserDataOperator::OPROW),
1081 cTx(ctx), F(f) {}
1082
1083 VectorDouble nF;
1084
1085 /**
1086 * \brief Integrate boundary condition
1087 * @param side local index of entity
1088 * @param type type of entity
1089 * @param data data on entity
1090 * @return error code
1091 */
1092 MoFEMErrorCode doWork(int side, EntityType type,
1093 EntitiesFieldData::EntData &data) {
1095 if (data.getFieldData().size() == 0)
1097 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1098 nF.resize(data.getIndices().size());
1099 nF.clear();
1100 int nb_gauss_pts = data.getN().size1();
1101 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1102 const double x = getCoordsAtGaussPts()(gg, 0);
1103 const double y = getCoordsAtGaussPts()(gg, 1);
1104 const double z = getCoordsAtGaussPts()(gg, 2);
1105
1106 double value;
1107 CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
1108 ;
1109 double w = getGaussPts()(2, gg) * 0.5;
1110 if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1111 noalias(nF) +=
1112 w * prod(data.getVectorN<3>(gg), getNormalsAtGaussPts(gg)) *
1113 value;
1114 } else {
1115 noalias(nF) += w * prod(data.getVectorN<3>(gg), getNormal()) * value;
1116 }
1117 }
1118 Vec f = (F != PETSC_NULL) ? F : getFEMethod()->ts_F;
1119 CHKERR VecSetValues(f, data.getIndices().size(), &data.getIndices()[0],
1121
1123 }
1124 };
1125
1126 /**
1127 * \brief Evaluate boundary conditions on fluxes.
1128 *
1129 * Note that Neumann boundary conditions here are essential. So it is opposite
1130 * what you find in displacement finite element method.
1131 *
1132
1133 * Here we have to solve for degrees of freedom on boundary such base
1134 functions
1135 * approximate flux.
1136
1137 *
1138 */
1142 OpEvaluateBcOnFluxes(MixTransportElement &ctx, const std::string flux_name)
1144 flux_name, UserDataOperator::OPROW),
1145 cTx(ctx) {}
1146
1147 MatrixDouble NN;
1148 VectorDouble Nf;
1150
1151 MoFEMErrorCode doWork(int side, EntityType type,
1152 EntitiesFieldData::EntData &data) {
1154 if (data.getFieldData().size() == 0)
1156 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1157 int nb_dofs = data.getFieldData().size();
1158 int nb_gauss_pts = data.getN().size1();
1159 if (3 * nb_dofs != static_cast<int>(data.getN().size2())) {
1160 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1161 "wrong number of dofs");
1162 }
1163 NN.resize(nb_dofs, nb_dofs);
1164 Nf.resize(nb_dofs);
1165 NN.clear();
1166 Nf.clear();
1167
1168 // Get normal vector. Note that when higher order geometry is set, then
1169 // face element could be curved, i.e. normal can be different at each
1170 // integration point.
1171 double *normal_ptr;
1172 if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1173 // HO geometry
1174 normal_ptr = &getNormalsAtGaussPts(0)[0];
1175 } else {
1176 // Linear geometry, i.e. constant normal on face
1177 normal_ptr = &getNormal()[0];
1178 }
1179 // set tensor from pointer
1180 FTensor::Tensor1<const double *, 3> t_normal(normal_ptr, &normal_ptr[1],
1181 &normal_ptr[2], 3);
1182 // get base functions
1183 auto t_n_hdiv_row = data.getFTensor1N<3>();
1184
1185 double nrm2 = 0;
1186
1187 // loop over integration points
1188 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1189
1190 // get integration point coordinates
1191 const double x = getCoordsAtGaussPts()(gg, 0);
1192 const double y = getCoordsAtGaussPts()(gg, 1);
1193 const double z = getCoordsAtGaussPts()(gg, 2);
1194
1195 // get flux on fece for given element handle and coordinates
1196 double flux;
1197 CHKERR cTx.getBcOnFluxes(fe_ent, x, y, z, flux);
1198 ;
1199 // get weight for integration rule
1200 double w = getGaussPts()(2, gg);
1201 if (gg == 0) {
1202 nrm2 = sqrt(t_normal(i) * t_normal(i));
1203 }
1204
1205 // set tensor of rank 0 to matrix NN elements
1206 // loop over base functions on rows and columns
1207 for (int ll = 0; ll != nb_dofs; ll++) {
1208 // get column on shape functions
1210 &data.getVectorN<3>(gg)(0, HVEC0),
1211 &data.getVectorN<3>(gg)(0, HVEC1),
1212 &data.getVectorN<3>(gg)(0, HVEC2), 3);
1213 for (int kk = 0; kk <= ll; kk++) {
1214 NN(ll, kk) += w * t_n_hdiv_row(i) * t_n_hdiv_col(i);
1215 ++t_n_hdiv_col;
1216 }
1217 // right hand side
1218 Nf[ll] += w * t_n_hdiv_row(i) * t_normal(i) * flux / nrm2;
1219 ++t_n_hdiv_row;
1220 }
1221
1222 // If HO geometry increment t_normal to next integration point
1223 if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1224 ++t_normal;
1225 nrm2 = sqrt(t_normal(i) * t_normal(i));
1226 }
1227 }
1228 // get global dofs indices on element
1229 cTx.bcIndices.insert(data.getIndices().begin(), data.getIndices().end());
1230 // factor matrix
1232 // solve local problem
1233 cholesky_solve(NN, Nf, ublas::lower());
1234
1235 // set solution to vector
1236 CHKERR VecSetOption(cTx.D0, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1237 CHKERR VecSetValues(cTx.D0, nb_dofs, &*data.getIndices().begin(),
1238 &*Nf.begin(), INSERT_VALUES);
1239 for (int dd = 0; dd != nb_dofs; ++dd)
1240 data.getFieldDofs()[dd]->getFieldData() = Nf[dd];
1241
1243 }
1244 };
1245
1246 /**
1247 * \brief Calculate values at integration points
1248 */
1251
1253
1255 const std::string val_name = "VALUES")
1256 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1257 val_name, UserDataOperator::OPROW),
1258 cTx(ctx) {}
1259
1261
1262 MoFEMErrorCode doWork(int side, EntityType type,
1263 EntitiesFieldData::EntData &data) {
1265 if (data.getFieldData().size() == 0)
1267
1268 int nb_gauss_pts = data.getN().size1();
1269 cTx.valuesAtGaussPts.resize(nb_gauss_pts);
1270 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1271 cTx.valuesAtGaussPts[gg] =
1272 inner_prod(trans(data.getN(gg)), data.getFieldData());
1273 }
1274
1276 }
1277 };
1278
1279 /**
1280 * \brief Calculate gradients of values at integration points
1281 */
1284
1286
1288 const std::string val_name = "VALUES")
1289 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1290 val_name, UserDataOperator::OPROW),
1291 cTx(ctx) {}
1293
1294 MoFEMErrorCode doWork(int side, EntityType type,
1295 EntitiesFieldData::EntData &data) {
1297 if (data.getFieldData().size() == 0)
1299 int nb_gauss_pts = data.getDiffN().size1();
1301 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1305 prod(trans(data.getDiffN(gg)), data.getFieldData());
1306 }
1308 }
1309 };
1310
1311 /**
1312 * \brief calculate flux at integration point
1313 */
1316
1318
1320 const std::string field_name)
1321 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1323 cTx(ctx) {}
1324
1325 VectorDouble divVec;
1326 MoFEMErrorCode doWork(int side, EntityType type,
1327 EntitiesFieldData::EntData &data) {
1329
1330 if (data.getFieldData().size() == 0)
1332 int nb_gauss_pts = data.getDiffN().size1();
1333 int nb_dofs = data.getFieldData().size();
1334 cTx.fluxesAtGaussPts.resize(3, nb_gauss_pts);
1335 cTx.divergenceAtGaussPts.resize(nb_gauss_pts);
1336 if (type == MBTRI && side == 0) {
1337 cTx.divergenceAtGaussPts.clear();
1338 cTx.fluxesAtGaussPts.clear();
1339 }
1340 divVec.resize(nb_dofs);
1341
1343 auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
1344
1345 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1346 for (auto &v : divVec) {
1347 v = t_base_diff_hdiv(i, i);
1348 ++t_base_diff_hdiv;
1349 }
1350 cTx.divergenceAtGaussPts[gg] += inner_prod(divVec, data.getFieldData());
1351 ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1352 cTx.fluxesAtGaussPts, gg);
1353 flux_at_gauss_pt +=
1354 prod(trans(data.getVectorN<3>(gg)), data.getFieldData());
1355 }
1357 }
1358 };
1359
1360 map<double, EntityHandle> errorMap;
1364
1365 /** \brief calculate error evaluator
1366 */
1367 struct OpError
1369
1371
1372 OpError(MixTransportElement &ctx, const std::string field_name)
1373 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1375 cTx(ctx) {}
1376 virtual ~OpError() {}
1377
1378 VectorDouble deltaFlux;
1379 MatrixDouble3by3 invK;
1380
1381 MoFEMErrorCode doWork(int side, EntityType type,
1382 EntitiesFieldData::EntData &data) {
1384 if (type != MBTET)
1386 invK.resize(3, 3, false);
1387 int nb_gauss_pts = data.getN().size1();
1388 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1389 double def_val = 0;
1390 Tag th_error_flux, th_error_div;
1391 CHKERR cTx.mField.get_moab().tag_get_handle(
1392 "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
1393 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1394 double *error_flux_ptr;
1395 CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1396 th_error_flux, &fe_ent, 1, (const void **)&error_flux_ptr);
1397
1398 CHKERR cTx.mField.get_moab().tag_get_handle(
1399 "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
1400 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1401 double *error_div_ptr;
1402 CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1403 th_error_div, &fe_ent, 1, (const void **)&error_div_ptr);
1404
1405 Tag th_error_jump;
1406 CHKERR cTx.mField.get_moab().tag_get_handle(
1407 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1408 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1409 double *error_jump_ptr;
1410 CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1411 th_error_jump, &fe_ent, 1, (const void **)&error_jump_ptr);
1412 *error_jump_ptr = 0;
1413
1414 /// characteristic size of the element
1415 const double h = pow(getVolume() * 12 / std::sqrt(2), (double)1 / 3);
1416
1417 for (int ff = 0; ff != 4; ff++) {
1418 EntityHandle face;
1419 CHKERR cTx.mField.get_moab().side_element(fe_ent, 2, ff, face);
1420 double *error_face_jump_ptr;
1421 CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1422 th_error_jump, &face, 1, (const void **)&error_face_jump_ptr);
1423 *error_face_jump_ptr = (1 / sqrt(h)) * sqrt(*error_face_jump_ptr);
1424 *error_face_jump_ptr = pow(*error_face_jump_ptr, 2);
1425 *error_jump_ptr += *error_face_jump_ptr;
1426 }
1427
1428 *error_flux_ptr = 0;
1429 *error_div_ptr = 0;
1430 deltaFlux.resize(3, false);
1431 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1432 double w = getGaussPts()(3, gg) * getVolume();
1433 const double x = getCoordsAtGaussPts()(gg, 0);
1434 const double y = getCoordsAtGaussPts()(gg, 1);
1435 const double z = getCoordsAtGaussPts()(gg, 2);
1436 double flux;
1437 CHKERR cTx.getSource(fe_ent, x, y, z, flux);
1438 ;
1439 *error_div_ptr += w * pow(cTx.divergenceAtGaussPts[gg] - flux, 2);
1440 CHKERR cTx.getResistivity(fe_ent, x, y, z, invK);
1441 ;
1442 ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1443 cTx.fluxesAtGaussPts, gg);
1446 noalias(deltaFlux) =
1448 *error_flux_ptr += w * inner_prod(deltaFlux, deltaFlux);
1449 }
1450 *error_div_ptr = h * sqrt(*error_div_ptr);
1451 *error_div_ptr = pow(*error_div_ptr, 2);
1452 cTx.errorMap[sqrt(*error_flux_ptr + *error_div_ptr + *error_jump_ptr)] =
1453 fe_ent;
1454 // Sum/Integrate all errors
1455 cTx.sumErrorFlux += *error_flux_ptr * getVolume();
1456 cTx.sumErrorDiv += *error_div_ptr * getVolume();
1457 // FIXME: Summation should be while skeleton is calculated
1458 cTx.sumErrorJump +=
1459 *error_jump_ptr * getVolume(); /// FIXME: this need to be fixed
1461 }
1462 };
1463
1464 /**
1465 * \brief calculate jump on entities
1466 */
1469
1470 /**
1471 * \brief volume element to get values from adjacent tets to face
1472 */
1473 VolumeElementForcesAndSourcesCoreOnSide volSideFe;
1474
1475 /** store values at integration point, key of the map is sense of face in
1476 * respect to adjacent tetrahedra
1477 */
1478 map<int, VectorDouble> valMap;
1479
1480 /**
1481 * \brief calculate values on adjacent tetrahedra to face
1482 */
1484 : public VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator {
1485 map<int, VectorDouble> &valMap;
1486 OpVolSide(map<int, VectorDouble> &val_map)
1487 : VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator(
1488 "VALUES", UserDataOperator::OPROW),
1489 valMap(val_map) {}
1490 MoFEMErrorCode doWork(int side, EntityType type,
1491 EntitiesFieldData::EntData &data) {
1493 if (data.getFieldData().size() == 0)
1495 int nb_gauss_pts = data.getN().size1();
1496 valMap[getFaceSense()].resize(nb_gauss_pts);
1497 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1498 valMap[getFaceSense()][gg] =
1499 inner_prod(trans(data.getN(gg)), data.getFieldData());
1500 }
1502 }
1503 };
1504
1506
1507 OpSkeleton(MixTransportElement &ctx, const std::string flux_name)
1509 flux_name, UserDataOperator::OPROW),
1510 volSideFe(ctx.mField), cTx(ctx) {
1511 volSideFe.getOpPtrVector().push_back(new OpSkeleton::OpVolSide(valMap));
1512 }
1513
1514 MoFEMErrorCode doWork(int side, EntityType type,
1515 EntitiesFieldData::EntData &data) {
1517
1518 if (type == MBTRI) {
1519 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1520
1521 double def_val = 0;
1522 Tag th_error_jump;
1523 CHKERR cTx.mField.get_moab().tag_get_handle(
1524 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1525 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1526 double *error_jump_ptr;
1527 CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1528 th_error_jump, &fe_ent, 1, (const void **)&error_jump_ptr);
1529 *error_jump_ptr = 0;
1530
1531 // check if this is essential boundary condition
1532 EntityHandle essential_bc_meshset =
1533 cTx.mField.get_finite_element_meshset("MIX_BCFLUX");
1534 if (cTx.mField.get_moab().contains_entities(essential_bc_meshset,
1535 &fe_ent, 1)) {
1536 // essential bc, np jump then, exit and go to next face
1538 }
1539
1540 // calculate values form adjacent tets
1541 valMap.clear();
1542 CHKERR loopSideVolumes("MIX", volSideFe);
1543 ;
1544
1545 int nb_gauss_pts = data.getN().size1();
1546
1547 // it is only one face, so it has to be bc natural boundary condition
1548 if (valMap.size() == 1) {
1549 if (static_cast<int>(valMap.begin()->second.size()) != nb_gauss_pts) {
1550 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1551 "wrong number of integration points");
1552 }
1553 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1554 const double x = getCoordsAtGaussPts()(gg, 0);
1555 const double y = getCoordsAtGaussPts()(gg, 1);
1556 const double z = getCoordsAtGaussPts()(gg, 2);
1557 double value;
1558 CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
1559 double w = getGaussPts()(2, gg);
1560 if (static_cast<int>(getNormalsAtGaussPts().size1()) ==
1561 nb_gauss_pts) {
1562 w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1563 } else {
1564 w *= getArea();
1565 }
1566 *error_jump_ptr += w * pow(value - valMap.begin()->second[gg], 2);
1567 }
1568 } else if (valMap.size() == 2) {
1569 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1570 double w = getGaussPts()(2, gg);
1571 if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1572 w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1573 } else {
1574 w *= getArea();
1575 }
1576 double delta = valMap.at(1)[gg] - valMap.at(-1)[gg];
1577 *error_jump_ptr += w * pow(delta, 2);
1578 }
1579 } else {
1580 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1581 "data inconsistency, wrong number of neighbors "
1582 "valMap.size() = %d",
1583 valMap.size());
1584 }
1585 }
1586
1588 }
1589 };
1590};
1591
1592} // namespace MixTransport
1593
1594#endif //_MIX_TRANPORT_ELEMENT_HPP_
1595
1596/***************************************************************************/ /**
1597 * \defgroup mofem_mix_transport_elem Mix transport element
1598 * \ingroup user_modules
1599 ******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
static const double eps
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
constexpr int order
@ F
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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 add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
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_row(const std::string &fe_name, const std::string name_row)=0
set field row 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
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 add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
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
FTensor::Index< 'i', SPACE_DIM > i
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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)
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
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)
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)
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
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