v0.6.9
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 [5] [6] [16] 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

Definition at line 78 of file MixTransportElement.hpp.

78  :
79  mField(m_field),
80  feVol(m_field),
81  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 208 of file MixTransportElement.hpp.

208  {
209 
211  //Fields
214 
215  //meshset consisting all entities in mesh
216  EntityHandle root_set = mField.get_moab().get_root_set();
217  //add entities to field
218  ierr = mField.add_ents_to_field_by_type(root_set,MBTET,fluxes); CHKERRG(ierr);
219  ierr = mField.add_ents_to_field_by_type(root_set,MBTET,values); CHKERRG(ierr);
220  ierr = mField.set_field_order(root_set,MBTET,fluxes,order+1); CHKERRG(ierr);
221  ierr = mField.set_field_order(root_set,MBTRI,fluxes,order+1); CHKERRG(ierr);
222  ierr = mField.set_field_order(root_set,MBTET,values,order); CHKERRG(ierr);
223 
224 
226  }
field with continuous normal traction
Definition: definitions.h:176
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
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=-1)=0
Add field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=-1)=0
Add entities to field meshsetThe lower dimension entities are added depending on the space type...
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=-1)=0
Set order approximation of the entities in the field.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:140
Construction of base is by Demkowicz .
Definition: definitions.h:143
field with C-1 continuity
Definition: definitions.h:178

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

233  {
235 
236 
237 
238 
239  // Set up volume element operators. Operators are used to calculate components
240  // of stiffness matrix & right hand side, in essence are used to do volume integrals over
241  // tetrahedral in this case.
242 
243  // Define element "MIX". Note that this element will work with fluxes_name and
244  // values_name. This reflect bilinear form for the problem
252 
253  // Define finite element to integrate over skeleton, we need that to evaluate error
254  ierr = mField.add_finite_element("MIX_SKELETON",MF_ZERO); CHKERRG(ierr);
255  ierr = mField.modify_finite_element_add_field_row("MIX_SKELETON",fluxes_name); CHKERRG(ierr);
256  ierr = mField.modify_finite_element_add_field_col("MIX_SKELETON",fluxes_name); CHKERRG(ierr);
257  ierr = mField.modify_finite_element_add_field_data("MIX_SKELETON",fluxes_name); CHKERRG(ierr);
258 
259  // In some cases you like to use HO geometry to describe shape of the bode, curved edges and faces, for
260  // example body is a sphere. HO geometry is approximated by a field, which can be hierarchical, so shape of
261  // the edges could be given by polynomial of arbitrary order.
262  //
263  // Check if field "mesh_nodals_positions" is defined, and if it is add that field to data of finite
264  // element. MoFEM will use that that to calculate Jacobian as result that geometry in nonlinear.
265  if(mField.check_field(mesh_nodals_positions)) {
266  ierr = mField.modify_finite_element_add_field_data("MIX",mesh_nodals_positions); CHKERRG(ierr);
267  ierr = mField.modify_finite_element_add_field_data("MIX_SKELETON",mesh_nodals_positions); CHKERRG(ierr);
268  }
269  // Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those BLOCKSETS
270  // and add them to "MIX" finite element. In addition get data form that meshset
271  // and set conductivity which is used to calculate fluxes from gradients of concentration
272  // or gradient of temperature, depending how you interpret variables.
274 
275  // cerr << *it << endl;
276 
277  Mat_Thermal temp_data;
278  ierr = it->getAttributeDataStructure(temp_data); CHKERRG(ierr);
279  setOfBlocks[it->getMeshsetId()].cOnductivity = temp_data.data.Conductivity;
280  setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
281  rval = mField.get_moab().get_entities_by_type(
282  it->meshset,MBTET,setOfBlocks[it->getMeshsetId()].tEts,true
283  ); CHKERRG(rval);
285  setOfBlocks[it->getMeshsetId()].tEts,MBTET,"MIX"
286  ); CHKERRG(ierr);
287 
288  Range skeleton;
289  rval = mField.get_moab().get_adjacencies(
290  setOfBlocks[it->getMeshsetId()].tEts,2,false,skeleton,moab::Interface::UNION
291  );
293  skeleton,MBTRI,"MIX_SKELETON"
294  ); CHKERRG(ierr);
295 
296  }
297 
298  // Define element to integrate natural boundary conditions, i.e. set values.
299  ierr = mField.add_finite_element("MIX_BCVALUE",MF_ZERO); CHKERRG(ierr);
300  ierr = mField.modify_finite_element_add_field_row("MIX_BCVALUE",fluxes_name); CHKERRG(ierr);
301  ierr = mField.modify_finite_element_add_field_col("MIX_BCVALUE",fluxes_name); CHKERRG(ierr);
302  ierr = mField.modify_finite_element_add_field_data("MIX_BCVALUE",fluxes_name); CHKERRG(ierr);
303  ierr = mField.modify_finite_element_add_field_data("MIX_BCVALUE",values_name); CHKERRG(ierr);
304  if(mField.check_field(mesh_nodals_positions)) {
305  ierr = mField.modify_finite_element_add_field_data("MIX_BCVALUE",mesh_nodals_positions); CHKERRG(ierr);
306  }
307 
308  // Define element to apply essential boundary conditions.
309  ierr = mField.add_finite_element("MIX_BCFLUX",MF_ZERO); CHKERRG(ierr);
310  ierr = mField.modify_finite_element_add_field_row("MIX_BCFLUX",fluxes_name); CHKERRG(ierr);
311  ierr = mField.modify_finite_element_add_field_col("MIX_BCFLUX",fluxes_name); CHKERRG(ierr);
312  ierr = mField.modify_finite_element_add_field_data("MIX_BCFLUX",fluxes_name); CHKERRG(ierr);
313  ierr = mField.modify_finite_element_add_field_data("MIX_BCFLUX",values_name); CHKERRG(ierr);
314  if(mField.check_field(mesh_nodals_positions)) {
315  ierr = mField.modify_finite_element_add_field_data("MIX_BCFLUX",mesh_nodals_positions); CHKERRG(ierr);
316  }
317 
319  }
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 MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
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
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
#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:236
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL)=0
add finite element
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

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

326  {
327 
329  //build field
331  // get tetrahedrons which has been build previously and now in so called garbage bit level
332  Range done_tets;
333  ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
334  BitRefLevel().set(0),BitRefLevel().set(),MBTET,done_tets
335  ); CHKERRG(ierr);
336  ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
337  BitRefLevel().set(BITREFLEVEL_SIZE-1),BitRefLevel().set(),MBTET,done_tets
338  ); CHKERRG(ierr);
339  // get tetrahedrons which belong to problem bit level
340  Range ref_tets;
341  ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
342  ref_level,BitRefLevel().set(),MBTET,ref_tets
343  ); CHKERRG(ierr);
344  ref_tets = subtract(ref_tets,done_tets);
345  ierr = mField.build_finite_elements("MIX",&ref_tets,2); CHKERRG(ierr);
346  // get triangles which has been build previously and now in so called garbage bit level
347  Range done_faces;
348  ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
349  BitRefLevel().set(0),BitRefLevel().set(),MBTRI,done_faces
350  ); CHKERRG(ierr);
351  ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
352  BitRefLevel().set(BITREFLEVEL_SIZE-1),BitRefLevel().set(),MBTRI,done_faces
353  ); CHKERRG(ierr);
354  // get triangles which belong to problem bit level
355  Range ref_faces;
356  ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
357  ref_level,BitRefLevel().set(),MBTRI,ref_faces
358  ); CHKERRG(ierr);
359  ref_faces = subtract(ref_faces,done_faces);
360  //build finite elements structures
361  ierr = mField.build_finite_elements("MIX_BCFLUX",&ref_faces,2); CHKERRG(ierr);
362  ierr = mField.build_finite_elements("MIX_BCVALUE",&ref_faces,2); CHKERRG(ierr);
363  ierr = mField.build_finite_elements("MIX_SKELETON",&ref_faces,2); CHKERRG(ierr);
364  //Build adjacencies of degrees of freedom and elements
365  ierr = mField.build_adjacencies(ref_level); CHKERRG(ierr);
366  //Define problem
368  //set refinement level for problem
370  // Add element to problem
372  ierr = mField.modify_problem_add_finite_element("MIX","MIX_SKELETON"); CHKERRG(ierr);
373  // Boundary conditions
374  ierr = mField.modify_problem_add_finite_element("MIX","MIX_BCFLUX"); CHKERRG(ierr);
375  ierr = mField.modify_problem_add_finite_element("MIX","MIX_BCVALUE"); CHKERRG(ierr);
376  //build problem
377 
378  ProblemsManager *prb_mng_ptr;
379  ierr = mField.getInterface(prb_mng_ptr); CHKERRG(ierr);
380  ierr = prb_mng_ptr->buildProblem("MIX",true); CHKERRG(ierr);
381  //mesh partitioning
382  //partition
383  ierr = prb_mng_ptr->partitionProblem("MIX"); CHKERRG(ierr);
384  ierr = prb_mng_ptr->partitionFiniteElements("MIX"); CHKERRG(ierr);
385  //what are ghost nodes, see Petsc Manual
386  ierr = prb_mng_ptr->partitionGhostDofs("MIX"); CHKERRG(ierr);
388  }
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=-1)=0
Add problem.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
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 CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
virtual MoFEMErrorCode build_finite_elements(int verb=-1)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
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:309
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problemif same finite element is solved using different level of refinements...
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
virtual MoFEMErrorCode build_fields(int verb=-1)=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=-1)=0
build adjacencies

◆ calculateResidual()

MoFEMErrorCode MixTransport::MixTransportElement::calculateResidual ( )

calculate residual

Definition at line 613 of file MixTransportElement.hpp.

613  {
614 
616  ierr = VecZeroEntries(F); CHKERRG(ierr);
617  ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
618  ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
619  ierr = VecAssemblyBegin(F); CHKERRG(ierr);
620  ierr = VecAssemblyEnd(F); CHKERRG(ierr);
621  //calculate residuals
622  feVol.getOpPtrVector().clear();
623  feVol.getOpPtrVector().push_back(new OpL2Source(*this,"VALUES",F));
624  feVol.getOpPtrVector().push_back(new OpFluxDivergenceAtGaussPts(*this,"FLUXES"));
625  feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this,"VALUES"));
626  feVol.getOpPtrVector().push_back(new OpDivTauU_HdivL2(*this,"FLUXES","VALUES",F));
627  feVol.getOpPtrVector().push_back(new OpTauDotSigma_HdivHdiv(*this,"FLUXES",PETSC_NULL,F));
628  feVol.getOpPtrVector().push_back(new OpVDivSigma_L2Hdiv(*this,"VALUES","FLUXES",PETSC_NULL,F));
630  feTri.getOpPtrVector().clear();
631  feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this,"FLUXES",F));
632  ierr = mField.loop_finite_elements("MIX","MIX_BCVALUE",feTri); CHKERRG(ierr);
633  ierr = VecGhostUpdateBegin(F,ADD_VALUES,SCATTER_REVERSE); CHKERRG(ierr);
634  ierr = VecGhostUpdateEnd(F,ADD_VALUES,SCATTER_REVERSE); CHKERRG(ierr);
635  ierr = VecAssemblyBegin(F); CHKERRG(ierr);
636  ierr = VecAssemblyEnd(F); CHKERRG(ierr);
637  // ierr = VecAXPY(F,-1.,D0); CHKERRG(ierr);
638  // ierr = MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULL,F); CHKERRG(ierr);
639  {
640  std::vector<int> ids;
641  ids.insert(ids.begin(),bcIndices.begin(),bcIndices.end());
642  std::vector<double> vals(ids.size(),0);
643  ierr = VecSetValues(F,ids.size(),&*ids.begin(),&*vals.begin(),INSERT_VALUES); CHKERRG(ierr);
644  ierr = VecAssemblyBegin(F); CHKERRG(ierr);
645  ierr = VecAssemblyEnd(F); CHKERRG(ierr);
646  }
647  {
648  double nrm2_F;
649  ierr = VecNorm(F,NORM_2,&nrm2_F); CHKERRG(ierr);
650  PetscPrintf(PETSC_COMM_WORLD,"nrm2_F = %6.4e\n",nrm2_F);
651  const double eps = 1e-8;
652  if(nrm2_F > eps) {
653  //SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"problem with residual");
654  }
655  }
657  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
MyVolumeFE feVol
Instance of volume element.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, MoFEMTypes bh=MF_EXIST, int verb=-1)=0
Make a loop over finite elements.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
MyTriFE feTri
Instance of surface element.
static const double eps

◆ createMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::createMatrices ( )

create matrices

Definition at line 484 of file MixTransportElement.hpp.

484  {
487  ierr = mField.getInterface<VecManager>()->vecCreateGhost("MIX",COL,&D); CHKERRG(ierr);
488  ierr = mField.getInterface<VecManager>()->vecCreateGhost("MIX",COL,&D0); CHKERRG(ierr);
489  ierr = mField.getInterface<VecManager>()->vecCreateGhost("MIX",ROW,&F); CHKERRG(ierr);
491  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
virtual MoFEMErrorCode MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=-1)=0
create Mat (MPIAIJ) for problem (collective)
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.

◆ destroyMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::destroyMatrices ( )

destroy matrices

Definition at line 692 of file MixTransportElement.hpp.

692  {
694  ierr = MatDestroy(&Aij); CHKERRG(ierr);
695  ierr = VecDestroy(&D); CHKERRG(ierr);
696  ierr = VecDestroy(&D0); CHKERRG(ierr);
697  ierr = VecDestroy(&F); CHKERRG(ierr);
699  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474

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

665  {
666 
668  errorMap.clear();
669  sumErrorFlux = 0;
670  sumErrorDiv = 0;
671  sumErrorJump = 0;
672  feTri.getOpPtrVector().clear();
673  feTri.getOpPtrVector().push_back(new OpSkeleton(*this,"FLUXES"));
674  ierr = mField.loop_finite_elements("MIX","MIX_SKELETON",feTri,0,mField.get_comm_size()); CHKERRG(ierr);
675  feVol.getOpPtrVector().clear();
676  feVol.getOpPtrVector().push_back(new OpFluxDivergenceAtGaussPts(*this,"FLUXES"));
677  feVol.getOpPtrVector().push_back(new OpValuesGradientAtGaussPts(*this,"VALUES"));
678  feVol.getOpPtrVector().push_back(new OpError(*this,"VALUES"));
680  const Problem *problem_ptr;
681  ierr = mField.get_problem("MIX",&problem_ptr); CHKERRG(ierr);
682  PetscPrintf(
683  mField.get_comm(),
684  "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error jump^2 = %6.4e error tot^2 = %6.4e\n",
685  problem_ptr->getNbDofsRow(),
687  );
689  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
virtual int get_comm_size() const =0
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
MyVolumeFE feVol
Instance of volume element.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, MoFEMTypes bh=MF_EXIST, int verb=-1)=0
Make a loop over finite elements.
map< double, EntityHandle > errorMap
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.
MyTriFE feTri
Instance of surface element.
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 182 of file MixTransportElement.hpp.

185  {
187  flux = 0;
189  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474

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

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

◆ getDirichletBCIndices()

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

get dof indices where essential boundary conditions are applied

Parameters
isindices
Returns
error code

Definition at line 100 of file MixTransportElement.hpp.

100  {
102  std::vector<int> ids;
103  ids.insert(ids.begin(),bcIndices.begin(),bcIndices.end());
104  IS is_local;
105  ierr = ISCreateGeneral(
106  mField.get_comm(),ids.size(),ids.empty()?PETSC_NULL:&ids[0],PETSC_COPY_VALUES,&is_local
107  ); CHKERRG(ierr);
108  ierr = ISAllGather(is_local,is); CHKERRG(ierr);
109  ierr = ISDestroy(&is_local); CHKERRG(ierr);
111  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
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.

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

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

◆ postProc()

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

Post process results.

Returns
error code

Definition at line 465 of file MixTransportElement.hpp.

465  {
466 
469  ierr = post_proc.generateReferenceElementMesh(); CHKERRG(ierr);
470  ierr = post_proc.addFieldValuesPostProc("VALUES"); CHKERRG(ierr);
471  ierr = post_proc.addFieldValuesGradientPostProc("VALUES"); CHKERRG(ierr);
472  ierr = post_proc.addFieldValuesPostProc("FLUXES"); CHKERRG(ierr);
473  // ierr = post_proc.addHdivFunctionsPostProc("FLUXES"); CHKERRG(ierr);
474  post_proc.getOpPtrVector().push_back(new OpPostProc(post_proc.postProcMesh,post_proc.mapGaussPts));
475  ierr = mField.loop_finite_elements("MIX","MIX",post_proc); CHKERRG(ierr);
476  ierr = post_proc.writeFile(out_file.c_str()); CHKERRG(ierr);
478  }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, MoFEMTypes bh=MF_EXIST, int verb=-1)=0
Make a loop over finite elements.
Post processing.

◆ solveLinearProblem()

MoFEMErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )

solve problem

Returns
error code

Definition at line 497 of file MixTransportElement.hpp.

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

Member Data Documentation

◆ Aij

Mat MixTransport::MixTransportElement::Aij

Definition at line 481 of file MixTransportElement.hpp.

◆ bcIndices

set<int> MixTransport::MixTransportElement::bcIndices

Definition at line 93 of file MixTransportElement.hpp.

◆ D

Vec MixTransport::MixTransportElement::D

Definition at line 480 of file MixTransportElement.hpp.

◆ D0

Vec MixTransport::MixTransportElement::D0

Definition at line 480 of file MixTransportElement.hpp.

◆ divergenceAtGaussPts

VectorDouble MixTransport::MixTransportElement::divergenceAtGaussPts

divergence at integration points on element

Definition at line 90 of file MixTransportElement.hpp.

◆ errorMap

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

Definition at line 1480 of file MixTransportElement.hpp.

◆ F

Vec MixTransport::MixTransportElement::F

Definition at line 480 of file MixTransportElement.hpp.

◆ feTri

MyTriFE MixTransport::MixTransportElement::feTri

Instance of surface element.

Definition at line 73 of file MixTransportElement.hpp.

◆ feVol

MyVolumeFE MixTransport::MixTransportElement::feVol

Instance of volume element.

Definition at line 60 of file MixTransportElement.hpp.

◆ fluxesAtGaussPts

MatrixDouble MixTransport::MixTransportElement::fluxesAtGaussPts

fluxes at integration points on element

Definition at line 91 of file MixTransportElement.hpp.

◆ mField

MoFEM::Interface& MixTransport::MixTransportElement::mField

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

◆ sumErrorDiv

double MixTransport::MixTransportElement::sumErrorDiv

Definition at line 1482 of file MixTransportElement.hpp.

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1481 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

Definition at line 1483 of file MixTransportElement.hpp.

◆ valuesAtGaussPts

VectorDouble MixTransport::MixTransportElement::valuesAtGaussPts

values at integration points on element

Definition at line 88 of file MixTransportElement.hpp.

◆ valuesGradientAtGaussPts

MatrixDouble MixTransport::MixTransportElement::valuesGradientAtGaussPts

gradients at integration points on element

Definition at line 89 of file MixTransportElement.hpp.


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