v0.14.0
Classes | Public Member Functions | Public Attributes | List of all members
MixTransport::MixTransportElement Struct Reference

Mix transport problem. More...

#include <users_modules/tutorials/cor-0to1/src/MixTransportElement.hpp>

Inheritance diagram for MixTransport::MixTransportElement:
[legend]
Collaboration diagram for MixTransport::MixTransportElement:
[legend]

Classes

struct  BlockData
 data for evaluation of het conductivity and heat capacity elements More...
 
struct  MyTriFE
 definition of surface element More...
 
struct  MyVolumeFE
 definition of volume element More...
 
struct  OpDivTauU_HdivL2
 Assemble \( \int_\mathcal{T} u \textrm{div}[\boldsymbol\tau] \textrm{d}\mathcal{T} \). More...
 
struct  OpError
 calculate error evaluator More...
 
struct  OpEvaluateBcOnFluxes
 Evaluate boundary conditions on fluxes. More...
 
struct  OpFluxDivergenceAtGaussPts
 calculate flux at integration point More...
 
struct  OpL2Source
 Calculate source therms, i.e. \(\int_\mathcal{T} f v \textrm{d}\mathcal{T}\). More...
 
struct  OpPostProc
 
struct  OpRhsBcOnValues
 calculate \( \int_\mathcal{S} {\boldsymbol\tau} \cdot \mathbf{n}u \textrm{d}\mathcal{S} \) More...
 
struct  OpSkeleton
 calculate jump on entities More...
 
struct  OpTauDotSigma_HdivHdiv
 Assemble \(\int_\mathcal{T} \mathbf{A} \boldsymbol\sigma \cdot \boldsymbol\tau \textrm{d}\mathcal{T}\). More...
 
struct  OpValuesAtGaussPts
 Calculate values at integration points. More...
 
struct  OpValuesGradientAtGaussPts
 Calculate gradients of values at integration points. More...
 
struct  OpVDivSigma_L2Hdiv
 \( \int_\mathcal{T} \textrm{div}[\boldsymbol\sigma] v \textrm{d}\mathcal{T} \) More...
 

Public Member Functions

 MixTransportElement (MoFEM::Interface &m_field)
 construction of this data structure More...
 
virtual ~MixTransportElement ()
 destructor More...
 
MoFEMErrorCode getDirichletBCIndices (IS *is)
 get dof indices where essential boundary conditions are applied More...
 
virtual MoFEMErrorCode getSource (const EntityHandle ent, const double x, const double y, const double z, double &flux)
 set source term More...
 
virtual MoFEMErrorCode getResistivity (const EntityHandle ent, const double x, const double y, const double z, MatrixDouble3by3 &inv_k)
 natural (Dirichlet) boundary conditions (set values) More...
 
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 More...
 
virtual MoFEMErrorCode getBcOnFluxes (const EntityHandle ent, const double x, const double y, const double z, double &flux)
 essential (Neumann) boundary condition (set fluxes) More...
 
MoFEMErrorCode addFields (const std::string &values, const std::string &fluxes, const int order)
 Add fields to database. More...
 
MoFEMErrorCode addFiniteElements (const std::string &fluxes_name, const std::string &values_name)
 add finite elements More...
 
MoFEMErrorCode buildProblem (BitRefLevel &ref_level)
 Build problem. More...
 
MoFEMErrorCode postProc (const string out_file)
 Post process results. More...
 
MoFEMErrorCode createMatrices ()
 create matrices More...
 
MoFEMErrorCode solveLinearProblem ()
 solve problem More...
 
MoFEMErrorCode calculateResidual ()
 calculate residual More...
 
MoFEMErrorCode evaluateError ()
 Calculate error on elements. More...
 
MoFEMErrorCode destroyMatrices ()
 destroy matrices More...
 

Public Attributes

MoFEM::InterfacemField
 
MyVolumeFE feVol
 Instance of volume element. More...
 
MyTriFE feTri
 Instance of surface element. More...
 
VectorDouble valuesAtGaussPts
 values at integration points on element More...
 
MatrixDouble valuesGradientAtGaussPts
 gradients at integration points on element More...
 
VectorDouble divergenceAtGaussPts
 divergence at integration points on element More...
 
MatrixDouble fluxesAtGaussPts
 fluxes at integration points on element More...
 
set< int > bcIndices
 
std::map< int, BlockDatasetOfBlocks
 maps block set id with appropriate BlockData More...
 
Vec D
 
Vec D0
 
Vec F
 
Mat Aij
 
map< double, EntityHandleerrorMap
 
double sumErrorFlux
 
double sumErrorDiv
 
double sumErrorJump
 

Detailed Description

Mix transport problem.

Note to solve this system you need to use direct solver or proper preconditioner for saddle problem.

It is based on [6] [7] [reviartthomas1996] 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

General problem have form,

\[ \mathbf{A} \boldsymbol\sigma + \textrm{grad}[u] = \mathbf{0} \; \textrm{on} \; \Omega \\ \textrm{div}[\boldsymbol\sigma] = f \; \textrm{on} \; \Omega \]

Definition at line 29 of file MixTransportElement.hpp.

Constructor & Destructor Documentation

◆ MixTransportElement()

MixTransport::MixTransportElement::MixTransportElement ( MoFEM::Interface m_field)
inline

construction of this data structure

Examples
UnsaturatedFlow.hpp.

Definition at line 65 of file MixTransportElement.hpp.

66  : mField(m_field), feVol(m_field), feTri(m_field){};

◆ ~MixTransportElement()

virtual MixTransport::MixTransportElement::~MixTransportElement ( )
inlinevirtual

destructor

Definition at line 71 of file MixTransportElement.hpp.

71 {}

Member Function Documentation

◆ addFields()

MoFEMErrorCode MixTransport::MixTransportElement::addFields ( const std::string &  values,
const std::string &  fluxes,
const int  order 
)
inline

Add fields to database.

Parameters
valuesname of the fields
fluxesname of filed for fluxes
orderorder of approximation
Returns
error code

Definition at line 189 of file MixTransportElement.hpp.

190  {
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
199  CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, fluxes);
200  CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, values);
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  }

◆ addFiniteElements()

MoFEMErrorCode MixTransport::MixTransportElement::addFiniteElements ( const std::string &  fluxes_name,
const std::string &  values_name 
)
inline

add finite elements

Definition at line 208 of file MixTransportElement.hpp.

209  {
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
228  CHKERR mField.add_finite_element("MIX_SKELETON", MF_ZERO);
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.
242  mField, BLOCKSET | MAT_THERMALSET, it)) {
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;
255  CHKERR mField.get_moab().get_adjacencies(
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.
263  CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
265  fluxes_name);
267  fluxes_name);
269  fluxes_name);
271  values_name);
272 
273  // Define element to apply essential boundary conditions.
274  CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
276  fluxes_name);
278  fluxes_name);
280  fluxes_name);
282  values_name);
284  }

◆ buildProblem()

MoFEMErrorCode MixTransport::MixTransportElement::buildProblem ( BitRefLevel &  ref_level)
inline

Build problem.

Parameters
ref_levelmesh refinement on which mesh problem you like to built.
Returns
error code

Definition at line 291 of file MixTransportElement.hpp.

291  {
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
327  CHKERR mField.build_adjacencies(ref_level);
328  // Define problem
330  // set refinement level for problem
332  // Add element to problem
334  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_SKELETON");
335  // Boundary conditions
336  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCFLUX");
337  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCVALUE");
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  }

◆ calculateResidual()

MoFEMErrorCode MixTransport::MixTransportElement::calculateResidual ( )
inline

calculate residual

Definition at line 612 of file MixTransportElement.hpp.

612  {
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();
621  CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
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);
636  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
637  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
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  }

◆ createMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::createMatrices ( )
inline

create matrices

Definition at line 474 of file MixTransportElement.hpp.

474  {
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  }

◆ destroyMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::destroyMatrices ( )
inline

destroy matrices

Definition at line 700 of file MixTransportElement.hpp.

700  {
702  CHKERR MatDestroy(&Aij);
703  ;
704  CHKERR VecDestroy(&D);
705  ;
706  CHKERR VecDestroy(&D0);
707  ;
708  CHKERR VecDestroy(&F);
709  ;
711  }

◆ evaluateError()

MoFEMErrorCode MixTransport::MixTransportElement::evaluateError ( )
inline

Calculate error on elements.

Todo:
this functions runs serial, in future should be parallel and work on distributed meshes.

Definition at line 670 of file MixTransportElement.hpp.

670  {
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,
680  feVol.getOpPtrVector().clear();
681  CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
682  feVol.getOpPtrVector().push_back(
683  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
684  feVol.getOpPtrVector().push_back(
685  new OpValuesGradientAtGaussPts(*this, "VALUES"));
686  feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
687  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
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,
697  }

◆ getBcOnFluxes()

virtual MoFEMErrorCode MixTransport::MixTransportElement::getBcOnFluxes ( const EntityHandle  ent,
const double  x,
const double  y,
const double  z,
double flux 
)
inlinevirtual

essential (Neumann) boundary condition (set fluxes)

Parameters
enthandle to finite element entity
xcoord
ycoord
zcoord
fluxreference to flux which is set by function
Returns
[description]

Reimplemented in MixTransport::UnsaturatedFlowElement, MyTransport, and MyTransport.

Definition at line 163 of file MixTransportElement.hpp.

165  {
167  flux = 0;
169  }

◆ getBcOnValues()

virtual MoFEMErrorCode MixTransport::MixTransportElement::getBcOnValues ( const EntityHandle  ent,
const int  gg,
const double  x,
const double  y,
const double  z,
double value 
)
inlinevirtual

evaluate natural (Dirichlet) boundary conditions

Parameters
ententity on which bc is evaluated
xcoordinate
ycoordinate
zcoordinate
valuevale
Returns
error code

Reimplemented in MixTransport::UnsaturatedFlowElement.

Definition at line 146 of file MixTransportElement.hpp.

148  {
150  value = 0;
152  }

◆ getDirichletBCIndices()

MoFEMErrorCode MixTransport::MixTransportElement::getDirichletBCIndices ( IS *  is)
inline

get dof indices where essential boundary conditions are applied

Parameters
isindices
Returns
error code

Definition at line 87 of file MixTransportElement.hpp.

87  {
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  }

◆ getResistivity()

virtual MoFEMErrorCode MixTransport::MixTransportElement::getResistivity ( const EntityHandle  ent,
const double  x,
const double  y,
const double  z,
MatrixDouble3by3 &  inv_k 
)
inlinevirtual

natural (Dirichlet) boundary conditions (set values)

Parameters
enthandle to finite element entity
xcoord
ycoord
zcoord
valuereference to value set by function
Returns
error code

Definition at line 126 of file MixTransportElement.hpp.

128  {
130  inv_k.clear();
131  for (int dd = 0; dd < 3; dd++) {
132  inv_k(dd, dd) = 1;
133  }
135  }

◆ getSource()

virtual MoFEMErrorCode MixTransport::MixTransportElement::getSource ( const EntityHandle  ent,
const double  x,
const double  y,
const double  z,
double flux 
)
inlinevirtual

set source term

Parameters
enthandle to entity on which function is evaluated
xcoord
ycoord
zcoord
fluxreference to source term set by function
Returns
error code

Reimplemented in MyTransport, and MyTransport.

Definition at line 109 of file MixTransportElement.hpp.

111  {
113  flux = 0;
115  }

◆ postProc()

MoFEMErrorCode MixTransport::MixTransportElement::postProc ( const string  out_file)
inline

Post process results.

Returns
error code

Definition at line 428 of file MixTransportElement.hpp.

428  {
430 
431  PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
432  mField);
433 
434  CHKERR AddHOOps<3, 3, 3>::add(post_proc.getOpPtrVector(), {HDIV, L2});
435 
436  auto values_ptr = boost::make_shared<VectorDouble>();
437  auto grad_ptr = boost::make_shared<MatrixDouble>();
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(
443  new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
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 
455  {{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
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  }

◆ solveLinearProblem()

MoFEMErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )
inline

solve problem

Returns
error code

Definition at line 488 of file MixTransportElement.hpp.

488  {
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();
511  CHKERR AddHOOps<2, 3, 3>::add(feTri.getOpPtrVector(), {HDIV});
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();
522  CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
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();
537  CHKERR AddHOOps<2, 3, 3>::add(feTri.getOpPtrVector(), {HDIV});
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);
544  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
545  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
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 
555  // CHKERR MatMultAdd(Aij,D0,F,F); ;
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  }

Member Data Documentation

◆ Aij

Mat MixTransport::MixTransportElement::Aij

Definition at line 471 of file MixTransportElement.hpp.

◆ bcIndices

set<int> MixTransport::MixTransportElement::bcIndices
Examples
UnsaturatedFlow.hpp.

Definition at line 80 of file MixTransportElement.hpp.

◆ D

Vec MixTransport::MixTransportElement::D

Definition at line 470 of file MixTransportElement.hpp.

◆ D0

Vec MixTransport::MixTransportElement::D0
Examples
UnsaturatedFlow.hpp.

Definition at line 470 of file MixTransportElement.hpp.

◆ divergenceAtGaussPts

VectorDouble MixTransport::MixTransportElement::divergenceAtGaussPts

divergence at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 77 of file MixTransportElement.hpp.

◆ errorMap

map<double, EntityHandle> MixTransport::MixTransportElement::errorMap

Definition at line 1360 of file MixTransportElement.hpp.

◆ F

Vec MixTransport::MixTransportElement::F

Definition at line 470 of file MixTransportElement.hpp.

◆ feTri

MyTriFE MixTransport::MixTransportElement::feTri

Instance of surface element.

Definition at line 60 of file MixTransportElement.hpp.

◆ feVol

MyVolumeFE MixTransport::MixTransportElement::feVol

Instance of volume element.

Definition at line 46 of file MixTransportElement.hpp.

◆ fluxesAtGaussPts

MatrixDouble MixTransport::MixTransportElement::fluxesAtGaussPts

fluxes at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 78 of file MixTransportElement.hpp.

◆ mField

MoFEM::Interface& MixTransport::MixTransportElement::mField
Examples
UnsaturatedFlow.hpp.

Definition at line 31 of file MixTransportElement.hpp.

◆ setOfBlocks

std::map<int, BlockData> MixTransport::MixTransportElement::setOfBlocks

maps block set id with appropriate BlockData

Definition at line 180 of file MixTransportElement.hpp.

◆ sumErrorDiv

double MixTransport::MixTransportElement::sumErrorDiv

Definition at line 1362 of file MixTransportElement.hpp.

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1361 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

Definition at line 1363 of file MixTransportElement.hpp.

◆ valuesAtGaussPts

VectorDouble MixTransport::MixTransportElement::valuesAtGaussPts

values at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 73 of file MixTransportElement.hpp.

◆ valuesGradientAtGaussPts

MatrixDouble MixTransport::MixTransportElement::valuesGradientAtGaussPts

gradients at integration points on element

Definition at line 75 of file MixTransportElement.hpp.


The documentation for this struct was generated from the following file:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::CoreInterface::loop_finite_elements
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.
OpError
Definition: initial_diffusion.cpp:61
MoFEM::CoreInterface::modify_problem_ref_level_set_bit
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
EntityHandle
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MixTransport::MixTransportElement::Aij
Mat Aij
Definition: MixTransportElement.hpp:471
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
MixTransport::MixTransportElement::sumErrorFlux
double sumErrorFlux
Definition: MixTransportElement.hpp:1361
MixTransport::MixTransportElement::getDirichletBCIndices
MoFEMErrorCode getDirichletBCIndices(IS *is)
get dof indices where essential boundary conditions are applied
Definition: MixTransportElement.hpp:87
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
order
constexpr int order
Definition: dg_projection.cpp:18
ROW
@ ROW
Definition: definitions.h:123
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MAT_THERMALSET
@ MAT_THERMALSET
block name is "MAT_THERMAL"
Definition: definitions.h:161
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::add_field
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.
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MixTransport::MixTransportElement::sumErrorDiv
double sumErrorDiv
Definition: MixTransportElement.hpp:1362
COL
@ COL
Definition: definitions.h:123
MixTransport::MixTransportElement::feTri
MyTriFE feTri
Instance of surface element.
Definition: MixTransportElement.hpp:60
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MixTransport::MixTransportElement::D
Vec D
Definition: MixTransportElement.hpp:470
MixTransport::MixTransportElement::F
Vec F
Definition: MixTransportElement.hpp:470
MixTransport::MixTransportElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: MixTransportElement.hpp:180
MixTransport::MixTransportElement::feVol
MyVolumeFE feVol
Instance of volume element.
Definition: MixTransportElement.hpp:46
MixTransport::MixTransportElement::D0
Vec D0
Definition: MixTransportElement.hpp:470
Range
FTensor::dd
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
MixTransport::MixTransportElement::bcIndices
set< int > bcIndices
Definition: MixTransportElement.hpp:80
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MixTransport::MixTransportElement::sumErrorJump
double sumErrorJump
Definition: MixTransportElement.hpp:1363
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
BITREFLEVEL_SIZE
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
MixTransport::MixTransportElement::errorMap
map< double, EntityHandle > errorMap
Definition: MixTransportElement.hpp:1360
MoFEM::CoreInterface::modify_problem_add_finite_element
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
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::CoreInterface::set_field_order
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.
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MixTransport::MixTransportElement::mField
MoFEM::Interface & mField
Definition: MixTransportElement.hpp:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698