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

Mix transport problemNote to solve this system you need to use direct solver or proper preconditioner for saddle problem. More...

#include <users_modules/basic_finite_elements/mix_transport/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, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
 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< intbcIndices
 
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 44 of file MixTransportElement.hpp.

Constructor & Destructor Documentation

◆ MixTransportElement()

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

construction of this data structure

Examples:
UnsaturatedFlow.hpp.

Definition at line 80 of file MixTransportElement.hpp.

81  : mField(m_field), feVol(m_field), feTri(m_field){};
MyVolumeFE feVol
Instance of volume element.
MyTriFE feTri
Instance of surface element.

◆ ~MixTransportElement()

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

destructor

Definition at line 86 of file MixTransportElement.hpp.

86 {}

Member Function Documentation

◆ addFields()

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

Add fields to database.

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

Definition at line 204 of file MixTransportElement.hpp.

205  {
207  // Fields
210 
211  // meshset consisting all entities in mesh
212  EntityHandle root_set = mField.get_moab().get_root_set();
213  // add entities to field
214  CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, fluxes);
215  CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, values);
216  CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
217  CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
218  CHKERR mField.set_field_order(root_set, MBTET, values, order);
220  }
field with continuous normal traction
Definition: definitions.h:170
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
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.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:142
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
field with C-1 continuity
Definition: definitions.h:171

◆ addFiniteElements()

MoFEMErrorCode MixTransport::MixTransportElement::addFiniteElements ( const std::string &  fluxes_name,
const std::string &  values_name,
const std::string  mesh_nodals_positions = "MESH_NODE_POSITIONS" 
)

add finite elements

Definition at line 223 of file MixTransportElement.hpp.

225  {
227 
228  // Set up volume element operators. Operators are used to calculate
229  // components of stiffness matrix & right hand side, in essence are used to
230  // do volume integrals over tetrahedral in this case.
231 
232  // Define element "MIX". Note that this element will work with fluxes_name
233  // and values_name. This reflect bilinear form for the problem
241 
242  // Define finite element to integrate over skeleton, we need that to
243  // evaluate error
244  CHKERR mField.add_finite_element("MIX_SKELETON", MF_ZERO);
246  fluxes_name);
248  fluxes_name);
250  fluxes_name);
251 
252  // In some cases you like to use HO geometry to describe shape of the bode,
253  // curved edges and faces, for example body is a sphere. HO geometry is
254  // approximated by a field, which can be hierarchical, so shape of the
255  // edges could be given by polynomial of arbitrary order.
256  //
257  // Check if field "mesh_nodals_positions" is defined, and if it is add that
258  // field to data of finite element. MoFEM will use that that to calculate
259  // Jacobian as result that geometry in nonlinear.
260  if (mField.check_field(mesh_nodals_positions)) {
262  mesh_nodals_positions);
264  mesh_nodals_positions);
265  }
266  // Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those
267  // BLOCKSETS and add them to "MIX" finite element. In addition get data form
268  // that meshset and set conductivity which is used to calculate fluxes from
269  // gradients of concentration or gradient of temperature, depending how you
270  // interpret variables.
272  mField, BLOCKSET | MAT_THERMALSET, it)) {
273 
274  // cerr << *it << endl;
275 
276  Mat_Thermal temp_data;
277  CHKERR it->getAttributeDataStructure(temp_data);
278  setOfBlocks[it->getMeshsetId()].cOnductivity =
279  temp_data.data.Conductivity;
280  setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
281  CHKERR mField.get_moab().get_entities_by_type(
282  it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
284  setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
285 
286  Range skeleton;
287  CHKERR mField.get_moab().get_adjacencies(
288  setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
289  moab::Interface::UNION);
291  "MIX_SKELETON");
292  }
293 
294  // Define element to integrate natural boundary conditions, i.e. set values.
295  CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
297  fluxes_name);
299  fluxes_name);
301  fluxes_name);
303  values_name);
304  if (mField.check_field(mesh_nodals_positions)) {
306  mesh_nodals_positions);
307  }
308 
309  // Define element to apply essential boundary conditions.
310  CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
312  fluxes_name);
314  fluxes_name);
316  fluxes_name);
318  values_name);
319  if (mField.check_field(mesh_nodals_positions)) {
321  mesh_nodals_positions);
322  }
323 
325  }
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
add entities to finite element
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
block name is "MAT_THERMAL"
Definition: definitions.h:221
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define CHKERR
Inline error check.
Definition: definitions.h:578
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 bool check_field(const std::string &name) const =0
check if field is in database
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
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403

◆ buildProblem()

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

Build problem.

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

Definition at line 332 of file MixTransportElement.hpp.

332  {
334  // build field
336  // get tetrahedrons which has been build previously and now in so called
337  // garbage bit level
338  Range done_tets;
339  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
340  BitRefLevel().set(0), BitRefLevel().set(), MBTET, done_tets);
341  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
342  BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTET,
343  done_tets);
344  // get tetrahedrons which belong to problem bit level
345  Range ref_tets;
346  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
347  ref_level, BitRefLevel().set(), MBTET, ref_tets);
348  ref_tets = subtract(ref_tets, done_tets);
349  CHKERR mField.build_finite_elements("MIX", &ref_tets, 2);
350  // get triangles which has been build previously and now in so called
351  // garbage bit level
352  Range done_faces;
353  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
354  BitRefLevel().set(0), BitRefLevel().set(), MBTRI, done_faces);
355  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
356  BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTRI,
357  done_faces);
358  // get triangles which belong to problem bit level
359  Range ref_faces;
360  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
361  ref_level, BitRefLevel().set(), MBTRI, ref_faces);
362  ref_faces = subtract(ref_faces, done_faces);
363  // build finite elements structures
364  CHKERR mField.build_finite_elements("MIX_BCFLUX", &ref_faces, 2);
365  CHKERR mField.build_finite_elements("MIX_BCVALUE", &ref_faces, 2);
366  CHKERR mField.build_finite_elements("MIX_SKELETON", &ref_faces, 2);
367  // Build adjacencies of degrees of freedom and elements
368  CHKERR mField.build_adjacencies(ref_level);
369  // Define problem
371  // set refinement level for problem
373  // Add element to problem
375  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_SKELETON");
376  // Boundary conditions
377  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCFLUX");
378  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCVALUE");
379  // build problem
380 
381  ProblemsManager *prb_mng_ptr;
382  CHKERR mField.getInterface(prb_mng_ptr);
383  CHKERR prb_mng_ptr->buildProblem("MIX", true);
384  // mesh partitioning
385  // partition
386  CHKERR prb_mng_ptr->partitionProblem("MIX");
387  CHKERR prb_mng_ptr->partitionFiniteElements("MIX");
388  // what are ghost nodes, see Petsc Manual
389  CHKERR prb_mng_ptr->partitionGhostDofs("MIX");
391  }
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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 ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:276
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
#define CHKERR
Inline error check.
Definition: definitions.h:578
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...

◆ calculateResidual()

MoFEMErrorCode MixTransport::MixTransportElement::calculateResidual ( )

calculate residual

Definition at line 625 of file MixTransportElement.hpp.

625  {
627  CHKERR VecZeroEntries(F);
628  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
629  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
630  CHKERR VecAssemblyBegin(F);
631  CHKERR VecAssemblyEnd(F);
632  // calculate residuals
633  feVol.getOpPtrVector().clear();
634  feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
635  feVol.getOpPtrVector().push_back(
636  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
637  feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
638  feVol.getOpPtrVector().push_back(
639  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
640  feVol.getOpPtrVector().push_back(
641  new OpTauDotSigma_HdivHdiv(*this, "FLUXES", PETSC_NULL, F));
642  feVol.getOpPtrVector().push_back(
643  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", PETSC_NULL, F));
644  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
645  feTri.getOpPtrVector().clear();
646  feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
647  CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
648  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
649  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
650  CHKERR VecAssemblyBegin(F);
651  CHKERR VecAssemblyEnd(F);
652  // CHKERR VecAXPY(F,-1.,D0);
653  // CHKERR MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULL,F);
654  {
655  std::vector<int> ids;
656  ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
657  std::vector<double> vals(ids.size(), 0);
658  CHKERR VecSetValues(F, ids.size(), &*ids.begin(), &*vals.begin(),
659  INSERT_VALUES);
660  CHKERR VecAssemblyBegin(F);
661  CHKERR VecAssemblyEnd(F);
662  }
663  {
664  double nrm2_F;
665  CHKERR VecNorm(F, NORM_2, &nrm2_F);
666  PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
667  const double eps = 1e-8;
668  if (nrm2_F > eps) {
669  // SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"problem with
670  // residual");
671  }
672  }
674  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
MyVolumeFE feVol
Instance of volume element.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, MoFEMTypes bh=MF_EXIST, int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:578
MyTriFE feTri
Instance of surface element.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
static const double eps

◆ createMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::createMatrices ( )

create matrices

Definition at line 488 of file MixTransportElement.hpp.

488  {
491  ;
492  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
493  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
494  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
496  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
#define CHKERR
Inline error check.
Definition: definitions.h:578
virtual MoFEMErrorCode MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=DEFAULT_VERBOSITY)=0
create Mat (MPIAIJ) for problem (collective)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403

◆ destroyMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::destroyMatrices ( )

destroy matrices

Definition at line 711 of file MixTransportElement.hpp.

711  {
713  CHKERR MatDestroy(&Aij);
714  ;
715  CHKERR VecDestroy(&D);
716  ;
717  CHKERR VecDestroy(&D0);
718  ;
719  CHKERR VecDestroy(&F);
720  ;
722  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403

◆ evaluateError()

MoFEMErrorCode MixTransport::MixTransportElement::evaluateError ( )

Calculate error on elements.

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

Definition at line 682 of file MixTransportElement.hpp.

682  {
684  errorMap.clear();
685  sumErrorFlux = 0;
686  sumErrorDiv = 0;
687  sumErrorJump = 0;
688  feTri.getOpPtrVector().clear();
689  feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
690  CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
692  feVol.getOpPtrVector().clear();
693  feVol.getOpPtrVector().push_back(
694  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
695  feVol.getOpPtrVector().push_back(
696  new OpValuesGradientAtGaussPts(*this, "VALUES"));
697  feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
698  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
700  const Problem *problem_ptr;
701  CHKERR mField.get_problem("MIX", &problem_ptr);
702  PetscPrintf(mField.get_comm(),
703  "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
704  "jump^2 = %6.4e error tot^2 = %6.4e\n",
705  problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
708  }
virtual int get_comm_size() const =0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
MyVolumeFE feVol
Instance of volume element.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, MoFEMTypes bh=MF_EXIST, int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
virtual MoFEMErrorCode get_problem(const std::string &problem_name, const Problem **problem_ptr) const =0
Get problem database (data structure)
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:578
MyTriFE feTri
Instance of surface element.
map< double, EntityHandle > errorMap
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
virtual MPI_Comm & get_comm() const =0

◆ getBcOnFluxes()

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

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

180  {
182  flux = 0;
184  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490

◆ getBcOnValues()

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

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

163  {
165  value = 0;
167  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490

◆ getDirichletBCIndices()

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

get dof indices where essential boundary conditions are applied

Parameters
isindices
Returns
error code

Definition at line 102 of file MixTransportElement.hpp.

102  {
104  std::vector<int> ids;
105  ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
106  IS is_local;
107  CHKERR ISCreateGeneral(mField.get_comm(), ids.size(),
108  ids.empty() ? PETSC_NULL : &ids[0],
109  PETSC_COPY_VALUES, &is_local);
110  CHKERR ISAllGather(is_local, is);
111  CHKERR ISDestroy(&is_local);
113  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
virtual MPI_Comm & get_comm() const =0

◆ getResistivity()

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

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

143  {
145  inv_k.clear();
146  for (int dd = 0; dd < 3; dd++) {
147  inv_k(dd, dd) = 1;
148  }
150  }
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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490

◆ getSource()

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

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

126  {
128  flux = 0;
130  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490

◆ postProc()

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

Post process results.

Returns
error code

Definition at line 469 of file MixTransportElement.hpp.

469  {
472  CHKERR post_proc.generateReferenceElementMesh();
473  CHKERR post_proc.addFieldValuesPostProc("VALUES");
474  CHKERR post_proc.addFieldValuesGradientPostProc("VALUES");
475  CHKERR post_proc.addFieldValuesPostProc("FLUXES");
476  // CHKERR post_proc.addHdivFunctionsPostProc("FLUXES");
477  post_proc.getOpPtrVector().push_back(
478  new OpPostProc(post_proc.postProcMesh, post_proc.mapGaussPts));
479  CHKERR mField.loop_finite_elements("MIX", "MIX", post_proc);
480  CHKERR post_proc.writeFile(out_file.c_str());
482  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, MoFEMTypes bh=MF_EXIST, int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
Post processing.
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403

◆ solveLinearProblem()

MoFEMErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )

solve problem

Returns
error code

Definition at line 502 of file MixTransportElement.hpp.

502  {
504  CHKERR MatZeroEntries(Aij);
505  CHKERR VecZeroEntries(F);
506  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
507  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
508  CHKERR VecZeroEntries(D0);
509  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
510  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
511  CHKERR VecZeroEntries(D);
512  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
513  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
514 
515  CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
516  "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
517 
518  // Calculate essential boundary conditions
519 
520  // clear essential bc indices, it could have dofs from other mesh refinement
521  bcIndices.clear();
522  // clear operator, just in case if some other operators are left on this
523  // element
524  feTri.getOpPtrVector().clear();
525  // set operator to calculate essential boundary conditions
526  feTri.getOpPtrVector().push_back(new OpEvaluateBcOnFluxes(*this, "FLUXES"));
527  CHKERR mField.loop_finite_elements("MIX", "MIX_BCFLUX", feTri);
528  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
529  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
530  CHKERR VecAssemblyBegin(D0);
531  CHKERR VecAssemblyEnd(D0);
532 
533  // set operators to calculate matrix and right hand side vectors
534  feVol.getOpPtrVector().clear();
535  feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
536  feVol.getOpPtrVector().push_back(
537  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
538  feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
539  feVol.getOpPtrVector().push_back(
540  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
541  feVol.getOpPtrVector().push_back(
542  new OpTauDotSigma_HdivHdiv(*this, "FLUXES", Aij, F));
543  feVol.getOpPtrVector().push_back(
544  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", Aij, F));
545  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
546  ;
547 
548  // calculate right hand side for natural boundary conditions
549  feTri.getOpPtrVector().clear();
550  feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
551  CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
552  ;
553 
554  // assemble matrices
555  CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
556  CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
557  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
558  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
559  CHKERR VecAssemblyBegin(F);
560  CHKERR VecAssemblyEnd(F);
561 
562  {
563  double nrm2_F;
564  CHKERR VecNorm(F, NORM_2, &nrm2_F);
565  PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
566  }
567 
568  // CHKERR MatMultAdd(Aij,D0,F,F); ;
569 
570  // for ksp solver vector is moved into rhs side
571  // for snes it is left ond the left
572  CHKERR VecScale(F, -1);
573 
574  IS essential_bc_ids;
575  CHKERR getDirichletBCIndices(&essential_bc_ids);
576  CHKERR MatZeroRowsColumnsIS(Aij, essential_bc_ids, 1, D0, F);
577  CHKERR ISDestroy(&essential_bc_ids);
578 
579  // {
580  // double norm;
581  // CHKERR MatNorm(Aij,NORM_FROBENIUS,&norm); ;
582  // PetscPrintf(PETSC_COMM_WORLD,"mat norm = %6.4e\n",norm);
583  // }
584 
585  {
586  double nrm2_F;
587  CHKERR VecNorm(F, NORM_2, &nrm2_F);
588  PetscPrintf(PETSC_COMM_WORLD, "With BC nrm2_F = %6.4e\n", nrm2_F);
589  }
590 
591  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
592  // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
593  // std::string wait;
594  // std::cin >> wait;
595 
596  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
597  // std::cin >> wait;
598 
599  // Solve
600  KSP solver;
601  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
602  CHKERR KSPSetOperators(solver, Aij, Aij);
603  CHKERR KSPSetFromOptions(solver);
604  CHKERR KSPSetUp(solver);
605  CHKERR KSPSolve(solver, F, D);
606  CHKERR KSPDestroy(&solver);
607 
608  {
609  double nrm2_D;
610  CHKERR VecNorm(D, NORM_2, &nrm2_D);
611  ;
612  PetscPrintf(PETSC_COMM_WORLD, "nrm2_D = %6.4e\n", nrm2_D);
613  }
614  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
615  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
616 
617  // copy data form vector on mesh
618  CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
619  "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
620 
622  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
MoFEMErrorCode getDirichletBCIndices(IS *is)
get dof indices where essential boundary conditions are applied
MyVolumeFE feVol
Instance of volume element.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, MoFEMTypes bh=MF_EXIST, int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:578
MyTriFE feTri
Instance of surface element.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403

Member Data Documentation

◆ Aij

Mat MixTransport::MixTransportElement::Aij
Examples:
UnsaturatedFlow.hpp.

Definition at line 485 of file MixTransportElement.hpp.

◆ bcIndices

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

Definition at line 95 of file MixTransportElement.hpp.

◆ D

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

Definition at line 484 of file MixTransportElement.hpp.

◆ D0

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

Definition at line 484 of file MixTransportElement.hpp.

◆ divergenceAtGaussPts

VectorDouble MixTransport::MixTransportElement::divergenceAtGaussPts

divergence at integration points on element

Examples:
UnsaturatedFlow.hpp.

Definition at line 92 of file MixTransportElement.hpp.

◆ errorMap

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

Definition at line 1395 of file MixTransportElement.hpp.

◆ F

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

Definition at line 484 of file MixTransportElement.hpp.

◆ feTri

MyTriFE MixTransport::MixTransportElement::feTri

Instance of surface element.

Definition at line 75 of file MixTransportElement.hpp.

◆ feVol

MyVolumeFE MixTransport::MixTransportElement::feVol

Instance of volume element.

Definition at line 61 of file MixTransportElement.hpp.

◆ fluxesAtGaussPts

MatrixDouble MixTransport::MixTransportElement::fluxesAtGaussPts

fluxes at integration points on element

Examples:
UnsaturatedFlow.hpp.

Definition at line 93 of file MixTransportElement.hpp.

◆ mField

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

Definition at line 46 of file MixTransportElement.hpp.

◆ setOfBlocks

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

maps block set id with appropriate BlockData

Definition at line 195 of file MixTransportElement.hpp.

◆ sumErrorDiv

double MixTransport::MixTransportElement::sumErrorDiv

Definition at line 1397 of file MixTransportElement.hpp.

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1396 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

Definition at line 1398 of file MixTransportElement.hpp.

◆ valuesAtGaussPts

VectorDouble MixTransport::MixTransportElement::valuesAtGaussPts

values at integration points on element

Examples:
UnsaturatedFlow.hpp.

Definition at line 88 of file MixTransportElement.hpp.

◆ valuesGradientAtGaussPts

MatrixDouble MixTransport::MixTransportElement::valuesGradientAtGaussPts

gradients at integration points on element

Definition at line 90 of file MixTransportElement.hpp.


The documentation for this struct was generated from the following file: