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

Mix transport problem. More...

#include <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::Interface & mField
 
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, BlockData > setOfBlocks
 maps block set id with appropriate BlockData More...
 
Vec D
 
Vec D0
 
Vec F
 
Mat Aij
 
map< double, EntityHandle > errorMap
 
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 34 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 70 of file MixTransportElement.hpp.

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

◆ ~MixTransportElement()

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

destructor

Definition at line 76 of file MixTransportElement.hpp.

76 {}

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 194 of file MixTransportElement.hpp.

195  {
196  MoFEMFunctionBegin;
197  // Fields
198  CHKERR mField.add_field(fluxes, HDIV, DEMKOWICZ_JACOBI_BASE, 1);
199  CHKERR mField.add_field(values, L2, AINSWORTH_LEGENDRE_BASE, 1);
200 
201  // meshset consisting all entities in mesh
202  EntityHandle root_set = mField.get_moab().get_root_set();
203  // add entities to field
204  CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, fluxes);
205  CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, values);
206  CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
207  CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
208  CHKERR mField.set_field_order(root_set, MBTET, values, order);
209  MoFEMFunctionReturn(0);
210  }

◆ addFiniteElements()

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

add finite elements

Definition at line 213 of file MixTransportElement.hpp.

214  {
215  MoFEMFunctionBegin;
216 
217  // Set up volume element operators. Operators are used to calculate
218  // components of stiffness matrix & right hand side, in essence are used to
219  // do volume integrals over tetrahedral in this case.
220 
221  // Define element "MIX". Note that this element will work with fluxes_name
222  // and values_name. This reflect bilinear form for the problem
223  CHKERR mField.add_finite_element("MIX", MF_ZERO);
224  CHKERR mField.modify_finite_element_add_field_row("MIX", fluxes_name);
225  CHKERR mField.modify_finite_element_add_field_col("MIX", fluxes_name);
226  CHKERR mField.modify_finite_element_add_field_row("MIX", values_name);
227  CHKERR mField.modify_finite_element_add_field_col("MIX", values_name);
228  CHKERR mField.modify_finite_element_add_field_data("MIX", fluxes_name);
229  CHKERR mField.modify_finite_element_add_field_data("MIX", values_name);
230 
231  // Define finite element to integrate over skeleton, we need that to
232  // evaluate error
233  CHKERR mField.add_finite_element("MIX_SKELETON", MF_ZERO);
234  CHKERR mField.modify_finite_element_add_field_row("MIX_SKELETON",
235  fluxes_name);
236  CHKERR mField.modify_finite_element_add_field_col("MIX_SKELETON",
237  fluxes_name);
238  CHKERR mField.modify_finite_element_add_field_data("MIX_SKELETON",
239  fluxes_name);
240 
241  // Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those
242  // BLOCKSETS and add them to "MIX" finite element. In addition get data form
243  // that meshset and set conductivity which is used to calculate fluxes from
244  // gradients of concentration or gradient of temperature, depending how you
245  // interpret variables.
246  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
247  mField, BLOCKSET | MAT_THERMALSET, it)) {
248 
249  Mat_Thermal temp_data;
250  CHKERR it->getAttributeDataStructure(temp_data);
251  setOfBlocks[it->getMeshsetId()].cOnductivity =
252  temp_data.data.Conductivity;
253  setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
254  CHKERR mField.get_moab().get_entities_by_type(
255  it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
256  CHKERR mField.add_ents_to_finite_element_by_type(
257  setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
258 
259  Range skeleton;
260  CHKERR mField.get_moab().get_adjacencies(
261  setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
262  moab::Interface::UNION);
263  CHKERR mField.add_ents_to_finite_element_by_type(skeleton, MBTRI,
264  "MIX_SKELETON");
265  }
266 
267  // Define element to integrate natural boundary conditions, i.e. set values.
268  CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
269  CHKERR mField.modify_finite_element_add_field_row("MIX_BCVALUE",
270  fluxes_name);
271  CHKERR mField.modify_finite_element_add_field_col("MIX_BCVALUE",
272  fluxes_name);
273  CHKERR mField.modify_finite_element_add_field_data("MIX_BCVALUE",
274  fluxes_name);
275  CHKERR mField.modify_finite_element_add_field_data("MIX_BCVALUE",
276  values_name);
277 
278  // Define element to apply essential boundary conditions.
279  CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
280  CHKERR mField.modify_finite_element_add_field_row("MIX_BCFLUX",
281  fluxes_name);
282  CHKERR mField.modify_finite_element_add_field_col("MIX_BCFLUX",
283  fluxes_name);
284  CHKERR mField.modify_finite_element_add_field_data("MIX_BCFLUX",
285  fluxes_name);
286  CHKERR mField.modify_finite_element_add_field_data("MIX_BCFLUX",
287  values_name);
288  MoFEMFunctionReturn(0);
289  }

◆ 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 296 of file MixTransportElement.hpp.

296  {
297  MoFEMFunctionBegin;
298  // build field
299  CHKERR mField.build_fields();
300  // get tetrahedrons which has been build previously and now in so called
301  // garbage bit level
302  Range done_tets;
303  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
304  BitRefLevel().set(0), BitRefLevel().set(), MBTET, done_tets);
305  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
306  BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTET,
307  done_tets);
308  // get tetrahedrons which belong to problem bit level
309  Range ref_tets;
310  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
311  ref_level, BitRefLevel().set(), MBTET, ref_tets);
312  ref_tets = subtract(ref_tets, done_tets);
313  CHKERR mField.build_finite_elements("MIX", &ref_tets, 2);
314  // get triangles which has been build previously and now in so called
315  // garbage bit level
316  Range done_faces;
317  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
318  BitRefLevel().set(0), BitRefLevel().set(), MBTRI, done_faces);
319  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
320  BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTRI,
321  done_faces);
322  // get triangles which belong to problem bit level
323  Range ref_faces;
324  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
325  ref_level, BitRefLevel().set(), MBTRI, ref_faces);
326  ref_faces = subtract(ref_faces, done_faces);
327  // build finite elements structures
328  CHKERR mField.build_finite_elements("MIX_BCFLUX", &ref_faces, 2);
329  CHKERR mField.build_finite_elements("MIX_BCVALUE", &ref_faces, 2);
330  CHKERR mField.build_finite_elements("MIX_SKELETON", &ref_faces, 2);
331  // Build adjacencies of degrees of freedom and elements
332  CHKERR mField.build_adjacencies(ref_level);
333  // Define problem
334  CHKERR mField.add_problem("MIX", MF_ZERO);
335  // set refinement level for problem
336  CHKERR mField.modify_problem_ref_level_set_bit("MIX", ref_level);
337  // Add element to problem
338  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX");
339  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_SKELETON");
340  // Boundary conditions
341  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCFLUX");
342  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCVALUE");
343  // build problem
344 
345  ProblemsManager *prb_mng_ptr;
346  CHKERR mField.getInterface(prb_mng_ptr);
347  CHKERR prb_mng_ptr->buildProblem("MIX", true);
348  // mesh partitioning
349  // partition
350  CHKERR prb_mng_ptr->partitionProblem("MIX");
351  CHKERR prb_mng_ptr->partitionFiniteElements("MIX");
352  // what are ghost nodes, see Petsc Manual
353  CHKERR prb_mng_ptr->partitionGhostDofs("MIX");
354  MoFEMFunctionReturn(0);
355  }

◆ calculateResidual()

MoFEMErrorCode MixTransport::MixTransportElement::calculateResidual ( )
inline

calculate residual

Definition at line 617 of file MixTransportElement.hpp.

617  {
618  MoFEMFunctionBegin;
619  CHKERR VecZeroEntries(F);
620  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
621  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
622  CHKERR VecAssemblyBegin(F);
623  CHKERR VecAssemblyEnd(F);
624  // calculate residuals
625  feVol.getOpPtrVector().clear();
626  CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
627  feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
628  feVol.getOpPtrVector().push_back(
629  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
630  feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
631  feVol.getOpPtrVector().push_back(
632  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
633  feVol.getOpPtrVector().push_back(
634  new OpTauDotSigma_HdivHdiv(*this, "FLUXES", PETSC_NULL, F));
635  feVol.getOpPtrVector().push_back(
636  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", PETSC_NULL, F));
637  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
638  feTri.getOpPtrVector().clear();
639  feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
640  CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
641  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
642  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
643  CHKERR VecAssemblyBegin(F);
644  CHKERR VecAssemblyEnd(F);
645  // CHKERR VecAXPY(F,-1.,D0);
646  // CHKERR MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULL,F);
647  {
648  std::vector<int> ids;
649  ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
650  std::vector<double> vals(ids.size(), 0);
651  CHKERR VecSetValues(F, ids.size(), &*ids.begin(), &*vals.begin(),
652  INSERT_VALUES);
653  CHKERR VecAssemblyBegin(F);
654  CHKERR VecAssemblyEnd(F);
655  }
656  {
657  double nrm2_F;
658  CHKERR VecNorm(F, NORM_2, &nrm2_F);
659  PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
660  const double eps = 1e-8;
661  if (nrm2_F > eps) {
662  // SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"problem with
663  // residual");
664  }
665  }
666  MoFEMFunctionReturn(0);
667  }

◆ createMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::createMatrices ( )
inline

create matrices

Definition at line 479 of file MixTransportElement.hpp.

479  {
480  MoFEMFunctionBegin;
481  CHKERR mField.getInterface<MatrixManager>()
482  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("MIX", &Aij);
483  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
484  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
485  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
486  MoFEMFunctionReturn(0);
487  }

◆ destroyMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::destroyMatrices ( )
inline

destroy matrices

Definition at line 705 of file MixTransportElement.hpp.

705  {
706  MoFEMFunctionBegin;
707  CHKERR MatDestroy(&Aij);
708  ;
709  CHKERR VecDestroy(&D);
710  ;
711  CHKERR VecDestroy(&D0);
712  ;
713  CHKERR VecDestroy(&F);
714  ;
715  MoFEMFunctionReturn(0);
716  }

◆ 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 675 of file MixTransportElement.hpp.

675  {
676  MoFEMFunctionBegin;
677  errorMap.clear();
678  sumErrorFlux = 0;
679  sumErrorDiv = 0;
680  sumErrorJump = 0;
681  feTri.getOpPtrVector().clear();
682  feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
683  CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
684  mField.get_comm_size());
685  feVol.getOpPtrVector().clear();
686  CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
687  feVol.getOpPtrVector().push_back(
688  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
689  feVol.getOpPtrVector().push_back(
690  new OpValuesGradientAtGaussPts(*this, "VALUES"));
691  feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
692  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
693  mField.get_comm_size());
694  const Problem *problem_ptr;
695  CHKERR mField.get_problem("MIX", &problem_ptr);
696  PetscPrintf(mField.get_comm(),
697  "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
698  "jump^2 = %6.4e error tot^2 = %6.4e\n",
699  problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
700  sumErrorJump, sumErrorFlux + sumErrorDiv + sumErrorJump);
701  MoFEMFunctionReturn(0);
702  }

◆ 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 168 of file MixTransportElement.hpp.

170  {
171  MoFEMFunctionBeginHot;
172  flux = 0;
173  MoFEMFunctionReturnHot(0);
174  }

◆ 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 151 of file MixTransportElement.hpp.

153  {
154  MoFEMFunctionBeginHot;
155  value = 0;
156  MoFEMFunctionReturnHot(0);
157  }

◆ 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 92 of file MixTransportElement.hpp.

92  {
93  MoFEMFunctionBegin;
94  std::vector<int> ids;
95  ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
96  IS is_local;
97  CHKERR ISCreateGeneral(mField.get_comm(), ids.size(),
98  ids.empty() ? PETSC_NULL : &ids[0],
99  PETSC_COPY_VALUES, &is_local);
100  CHKERR ISAllGather(is_local, is);
101  CHKERR ISDestroy(&is_local);
102  MoFEMFunctionReturn(0);
103  }

◆ 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 131 of file MixTransportElement.hpp.

133  {
134  MoFEMFunctionBeginHot;
135  inv_k.clear();
136  for (int dd = 0; dd < 3; dd++) {
137  inv_k(dd, dd) = 1;
138  }
139  MoFEMFunctionReturnHot(0);
140  }

◆ 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 114 of file MixTransportElement.hpp.

116  {
117  MoFEMFunctionBeginHot;
118  flux = 0;
119  MoFEMFunctionReturnHot(0);
120  }

◆ postProc()

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

Post process results.

Returns
error code

Definition at line 433 of file MixTransportElement.hpp.

433  {
434  MoFEMFunctionBegin;
435 
436  PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
437  mField);
438 
439  CHKERR AddHOOps<3, 3, 3>::add(post_proc.getOpPtrVector(), {HDIV, L2});
440 
441  auto values_ptr = boost::make_shared<VectorDouble>();
442  auto grad_ptr = boost::make_shared<MatrixDouble>();
443  auto flux_ptr = boost::make_shared<MatrixDouble>();
444 
445  post_proc.getOpPtrVector().push_back(
446  new OpCalculateScalarFieldValues("VALUES", values_ptr));
447  post_proc.getOpPtrVector().push_back(
448  new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
449  post_proc.getOpPtrVector().push_back(
450  new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
451 
452  using OpPPMap = OpPostProcMapInMoab<3, 3>;
453 
454  post_proc.getOpPtrVector().push_back(
455 
456  new OpPPMap(post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
457 
458  {{"VALUES", values_ptr}},
459 
460  {{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
461 
462  {}, {}
463 
464  ));
465 
466  post_proc.getOpPtrVector().push_back(new OpPostProc(
467  post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
468 
469  CHKERR mField.loop_finite_elements("MIX", "MIX", post_proc);
470 
471  CHKERR post_proc.writeFile(out_file.c_str());
472  MoFEMFunctionReturn(0);
473  }

◆ solveLinearProblem()

MoFEMErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )
inline

solve problem

Returns
error code

Definition at line 493 of file MixTransportElement.hpp.

493  {
494  MoFEMFunctionBegin;
495  CHKERR MatZeroEntries(Aij);
496  CHKERR VecZeroEntries(F);
497  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
498  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
499  CHKERR VecZeroEntries(D0);
500  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
501  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
502  CHKERR VecZeroEntries(D);
503  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
504  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
505 
506  CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
507  "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
508 
509  // Calculate essential boundary conditions
510 
511  // clear essential bc indices, it could have dofs from other mesh refinement
512  bcIndices.clear();
513  // clear operator, just in case if some other operators are left on this
514  // element
515  feTri.getOpPtrVector().clear();
516  CHKERR AddHOOps<2, 3, 3>::add(feTri.getOpPtrVector(), {HDIV});
517  // set operator to calculate essential boundary conditions
518  feTri.getOpPtrVector().push_back(new OpEvaluateBcOnFluxes(*this, "FLUXES"));
519  CHKERR mField.loop_finite_elements("MIX", "MIX_BCFLUX", feTri);
520  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
521  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
522  CHKERR VecAssemblyBegin(D0);
523  CHKERR VecAssemblyEnd(D0);
524 
525  // set operators to calculate matrix and right hand side vectors
526  feVol.getOpPtrVector().clear();
527  CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
528  feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
529  feVol.getOpPtrVector().push_back(
530  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
531  feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
532  feVol.getOpPtrVector().push_back(
533  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
534  feVol.getOpPtrVector().push_back(
535  new OpTauDotSigma_HdivHdiv(*this, "FLUXES", Aij, F));
536  feVol.getOpPtrVector().push_back(
537  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", Aij, F));
538  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
539 
540  // calculate right hand side for natural boundary conditions
541  feTri.getOpPtrVector().clear();
542  CHKERR AddHOOps<2, 3, 3>::add(feTri.getOpPtrVector(), {HDIV});
543  feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
544  CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
545 
546  // assemble matrices
547  CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
548  CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
549  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
550  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
551  CHKERR VecAssemblyBegin(F);
552  CHKERR VecAssemblyEnd(F);
553 
554  {
555  double nrm2_F;
556  CHKERR VecNorm(F, NORM_2, &nrm2_F);
557  PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
558  }
559 
560  // CHKERR MatMultAdd(Aij,D0,F,F); ;
561 
562  // for ksp solver vector is moved into rhs side
563  // for snes it is left ond the left
564  CHKERR VecScale(F, -1);
565 
566  IS essential_bc_ids;
567  CHKERR getDirichletBCIndices(&essential_bc_ids);
568  CHKERR MatZeroRowsColumnsIS(Aij, essential_bc_ids, 1, D0, F);
569  CHKERR ISDestroy(&essential_bc_ids);
570 
571  // {
572  // double norm;
573  // CHKERR MatNorm(Aij,NORM_FROBENIUS,&norm); ;
574  // PetscPrintf(PETSC_COMM_WORLD,"mat norm = %6.4e\n",norm);
575  // }
576 
577  {
578  double nrm2_F;
579  CHKERR VecNorm(F, NORM_2, &nrm2_F);
580  PetscPrintf(PETSC_COMM_WORLD, "With BC nrm2_F = %6.4e\n", nrm2_F);
581  }
582 
583  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
584  // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
585  // std::string wait;
586  // std::cin >> wait;
587 
588  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
589  // std::cin >> wait;
590 
591  // Solve
592  KSP solver;
593  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
594  CHKERR KSPSetOperators(solver, Aij, Aij);
595  CHKERR KSPSetFromOptions(solver);
596  CHKERR KSPSetUp(solver);
597  CHKERR KSPSolve(solver, F, D);
598  CHKERR KSPDestroy(&solver);
599 
600  {
601  double nrm2_D;
602  CHKERR VecNorm(D, NORM_2, &nrm2_D);
603  ;
604  PetscPrintf(PETSC_COMM_WORLD, "nrm2_D = %6.4e\n", nrm2_D);
605  }
606  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
607  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
608 
609  // copy data form vector on mesh
610  CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
611  "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
612 
613  MoFEMFunctionReturn(0);
614  }

Member Data Documentation

◆ Aij

Mat MixTransport::MixTransportElement::Aij

Definition at line 476 of file MixTransportElement.hpp.

◆ bcIndices

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

Definition at line 85 of file MixTransportElement.hpp.

◆ D

Vec MixTransport::MixTransportElement::D

Definition at line 475 of file MixTransportElement.hpp.

◆ D0

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

Definition at line 475 of file MixTransportElement.hpp.

◆ divergenceAtGaussPts

VectorDouble MixTransport::MixTransportElement::divergenceAtGaussPts

divergence at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 82 of file MixTransportElement.hpp.

◆ errorMap

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

Definition at line 1365 of file MixTransportElement.hpp.

◆ F

Vec MixTransport::MixTransportElement::F

Definition at line 475 of file MixTransportElement.hpp.

◆ feTri

MyTriFE MixTransport::MixTransportElement::feTri

Instance of surface element.

Definition at line 65 of file MixTransportElement.hpp.

◆ feVol

MyVolumeFE MixTransport::MixTransportElement::feVol

Instance of volume element.

Definition at line 51 of file MixTransportElement.hpp.

◆ fluxesAtGaussPts

MatrixDouble MixTransport::MixTransportElement::fluxesAtGaussPts

fluxes at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 83 of file MixTransportElement.hpp.

◆ mField

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

Definition at line 36 of file MixTransportElement.hpp.

◆ setOfBlocks

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

maps block set id with appropriate BlockData

Definition at line 185 of file MixTransportElement.hpp.

◆ sumErrorDiv

double MixTransport::MixTransportElement::sumErrorDiv

Definition at line 1367 of file MixTransportElement.hpp.

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1366 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

Definition at line 1368 of file MixTransportElement.hpp.

◆ valuesAtGaussPts

VectorDouble MixTransport::MixTransportElement::valuesAtGaussPts

values at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 78 of file MixTransportElement.hpp.

◆ valuesGradientAtGaussPts

MatrixDouble MixTransport::MixTransportElement::valuesGradientAtGaussPts

gradients at integration points on element

Definition at line 80 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:460
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.
MoFEM::ProblemsManager::buildProblem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
Definition: ProblemsManager.cpp:87
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::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
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:476
_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.
MoFEM::Mat_Thermal::data
_data_ data
Definition: MaterialBlocks.hpp:217
MixTransport::MixTransportElement::sumErrorFlux
double sumErrorFlux
Definition: MixTransportElement.hpp:1366
MixTransport::MixTransportElement::getDirichletBCIndices
MoFEMErrorCode getDirichletBCIndices(IS *is)
get dof indices where essential boundary conditions are applied
Definition: MixTransportElement.hpp:92
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
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
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
ROW
@ ROW
Definition: definitions.h:136
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:548
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:174
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition: UserDataOperators.hpp:2115
MoFEM::Problem::getNbDofsRow
DofIdx getNbDofsRow() const
Definition: ProblemsMultiIndices.hpp:376
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MixTransport::MixTransportElement::sumErrorDiv
double sumErrorDiv
Definition: MixTransportElement.hpp:1367
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
COL
@ COL
Definition: definitions.h:136
MixTransport::MixTransportElement::feTri
MyTriFE feTri
Instance of surface element.
Definition: MixTransportElement.hpp:65
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:670
MixTransport::MixTransportElement::D
Vec D
Definition: MixTransportElement.hpp:475
MixTransport::MixTransportElement::F
Vec F
Definition: MixTransportElement.hpp:475
MixTransport::MixTransportElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: MixTransportElement.hpp:185
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_field)=0
set finite element field data
MoFEM::ProblemsManager::partitionFiniteElements
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
Definition: ProblemsManager.cpp:2167
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
MixTransport::MixTransportElement::feVol
MyVolumeFE feVol
Instance of volume element.
Definition: MixTransportElement.hpp:51
MoFEM::MatrixManager
Matrix manager is used to build and partition problems.
Definition: MatrixManager.hpp:21
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MixTransport::MixTransportElement::D0
Vec D0
Definition: MixTransportElement.hpp:475
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:111
MixTransport::MixTransportElement::bcIndices
set< int > bcIndices
Definition: MixTransportElement.hpp:85
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MixTransport::MixTransportElement::sumErrorJump
double sumErrorJump
Definition: MixTransportElement.hpp:1368
MoFEM::Mat_Thermal
Thermal material data structure.
Definition: MaterialBlocks.hpp:201
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
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:232
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
MixTransport::MixTransportElement::errorMap
map< double, EntityHandle > errorMap
Definition: MixTransportElement.hpp:1365
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
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:453
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
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::ProblemsManager::partitionGhostDofs
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
Definition: ProblemsManager.cpp:2339
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:36
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
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.
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::ProblemsManager::partitionProblem
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
Definition: ProblemsManager.cpp:1683
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
Generated by Doxygen 1.8.17 and hosted at University of Glasgow