v0.9.1
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< 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 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:178
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:482
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:151
#define CHKERR
Inline error check.
Definition: definitions.h:601
constexpr int order
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
field with C-1 continuity
Definition: definitions.h:179

◆ 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  Mat_Thermal temp_data;
275  CHKERR it->getAttributeDataStructure(temp_data);
276  setOfBlocks[it->getMeshsetId()].cOnductivity =
277  temp_data.data.Conductivity;
278  setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
279  CHKERR mField.get_moab().get_entities_by_type(
280  it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
282  setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
283 
284  Range skeleton;
285  CHKERR mField.get_moab().get_adjacencies(
286  setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
287  moab::Interface::UNION);
289  "MIX_SKELETON");
290  }
291 
292  // Define element to integrate natural boundary conditions, i.e. set values.
293  CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
295  fluxes_name);
297  fluxes_name);
299  fluxes_name);
301  values_name);
302  if (mField.check_field(mesh_nodals_positions)) {
304  mesh_nodals_positions);
305  }
306 
307  // Define element to apply essential boundary conditions.
308  CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
310  fluxes_name);
312  fluxes_name);
314  fluxes_name);
316  values_name);
317  if (mField.check_field(mesh_nodals_positions)) {
319  mesh_nodals_positions);
320  }
321 
323  }
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:482
#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:229
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:601
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:412

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

330  {
332  // build field
334  // get tetrahedrons which has been build previously and now in so called
335  // garbage bit level
336  Range done_tets;
337  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
338  BitRefLevel().set(0), BitRefLevel().set(), MBTET, done_tets);
339  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
340  BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTET,
341  done_tets);
342  // get tetrahedrons which belong to problem bit level
343  Range ref_tets;
344  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
345  ref_level, BitRefLevel().set(), MBTET, ref_tets);
346  ref_tets = subtract(ref_tets, done_tets);
347  CHKERR mField.build_finite_elements("MIX", &ref_tets, 2);
348  // get triangles which has been build previously and now in so called
349  // garbage bit level
350  Range done_faces;
351  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
352  BitRefLevel().set(0), BitRefLevel().set(), MBTRI, done_faces);
353  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
354  BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTRI,
355  done_faces);
356  // get triangles which belong to problem bit level
357  Range ref_faces;
358  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
359  ref_level, BitRefLevel().set(), MBTRI, ref_faces);
360  ref_faces = subtract(ref_faces, done_faces);
361  // build finite elements structures
362  CHKERR mField.build_finite_elements("MIX_BCFLUX", &ref_faces, 2);
363  CHKERR mField.build_finite_elements("MIX_BCVALUE", &ref_faces, 2);
364  CHKERR mField.build_finite_elements("MIX_SKELETON", &ref_faces, 2);
365  // Build adjacencies of degrees of freedom and elements
366  CHKERR mField.build_adjacencies(ref_level);
367  // Define problem
369  // set refinement level for problem
371  // Add element to problem
373  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_SKELETON");
374  // Boundary conditions
375  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCFLUX");
376  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCVALUE");
377  // build problem
378 
379  ProblemsManager *prb_mng_ptr;
380  CHKERR mField.getInterface(prb_mng_ptr);
381  CHKERR prb_mng_ptr->buildProblem("MIX", true);
382  // mesh partitioning
383  // partition
384  CHKERR prb_mng_ptr->partitionProblem("MIX");
385  CHKERR prb_mng_ptr->partitionFiniteElements("MIX");
386  // what are ghost nodes, see Petsc Manual
387  CHKERR prb_mng_ptr->partitionGhostDofs("MIX");
389  }
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:482
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:284
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:601
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
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:412
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 623 of file MixTransportElement.hpp.

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

◆ createMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::createMatrices ( )

create matrices

Definition at line 486 of file MixTransportElement.hpp.

486  {
488  CHKERR mField.getInterface<MatrixManager>()
489  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("MIX", &Aij);
490  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
491  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
492  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
494  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:601
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ destroyMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::destroyMatrices ( )

destroy matrices

Definition at line 709 of file MixTransportElement.hpp.

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

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

680  {
682  errorMap.clear();
683  sumErrorFlux = 0;
684  sumErrorDiv = 0;
685  sumErrorJump = 0;
686  feTri.getOpPtrVector().clear();
687  feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
688  CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
690  feVol.getOpPtrVector().clear();
691  feVol.getOpPtrVector().push_back(
692  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
693  feVol.getOpPtrVector().push_back(
694  new OpValuesGradientAtGaussPts(*this, "VALUES"));
695  feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
696  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
698  const Problem *problem_ptr;
699  CHKERR mField.get_problem("MIX", &problem_ptr);
700  PetscPrintf(mField.get_comm(),
701  "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
702  "jump^2 = %6.4e error tot^2 = %6.4e\n",
703  problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
706  }
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:482
MyVolumeFE feVol
Instance of volume element.
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:601
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:412
virtual MPI_Comm & get_comm() const =0
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)

◆ 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:506
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513

◆ 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:506
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513

◆ 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:482
#define CHKERR
Inline error check.
Definition: definitions.h:601
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
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  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
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

◆ 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:506
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513

◆ postProc()

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

Post process results.

Returns
error code

Definition at line 467 of file MixTransportElement.hpp.

467  {
470  CHKERR post_proc.generateReferenceElementMesh();
471  CHKERR post_proc.addFieldValuesPostProc("VALUES");
472  CHKERR post_proc.addFieldValuesGradientPostProc("VALUES");
473  CHKERR post_proc.addFieldValuesPostProc("FLUXES");
474  // CHKERR post_proc.addHdivFunctionsPostProc("FLUXES");
475  post_proc.getOpPtrVector().push_back(
476  new OpPostProc(post_proc.postProcMesh, post_proc.mapGaussPts));
477  CHKERR mField.loop_finite_elements("MIX", "MIX", post_proc);
478  CHKERR post_proc.writeFile(out_file.c_str());
480  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
Post processing.
#define CHKERR
Inline error check.
Definition: definitions.h:601
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)

◆ solveLinearProblem()

MoFEMErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )

solve problem

Returns
error code

Definition at line 500 of file MixTransportElement.hpp.

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

Member Data Documentation

◆ Aij

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

Definition at line 483 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 482 of file MixTransportElement.hpp.

◆ D0

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

Definition at line 482 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 1392 of file MixTransportElement.hpp.

◆ F

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

Definition at line 482 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 1394 of file MixTransportElement.hpp.

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1393 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

Definition at line 1395 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: