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

Mix transport problem. More...

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

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

Classes

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

Public Member Functions

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

Public Attributes

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

Detailed Description

Mix transport problem.

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

It is based on [6] [7] [reviartthomas1996] https://www.researchgate.net/profile/Richard_Falk/publication/226454406_Differential_Complexes_and_Stability_of_Finite_Element_Methods_I._The_de_Rham_Complex/links/02e7e5214f0426ff77000000.pdf

General problem have form,

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

Definition at line 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){};
MyTriFE feTri
Instance of surface element.
MyVolumeFE feVol
Instance of volume 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  }
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
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.
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 moab::Interface & get_moab()=0

◆ addFiniteElements()

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

add finite elements

Definition at line 223 of file MixTransportElement.hpp.

224  {
226 
227  // Set up volume element operators. Operators are used to calculate
228  // components of stiffness matrix & right hand side, in essence are used to
229  // do volume integrals over tetrahedral in this case.
230 
231  // Define element "MIX". Note that this element will work with fluxes_name
232  // and values_name. This reflect bilinear form for the problem
240 
241  // Define finite element to integrate over skeleton, we need that to
242  // evaluate error
243  CHKERR mField.add_finite_element("MIX_SKELETON", MF_ZERO);
245  fluxes_name);
247  fluxes_name);
249  fluxes_name);
250 
251  // Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those
252  // BLOCKSETS and add them to "MIX" finite element. In addition get data form
253  // that meshset and set conductivity which is used to calculate fluxes from
254  // gradients of concentration or gradient of temperature, depending how you
255  // interpret variables.
257  mField, BLOCKSET | MAT_THERMALSET, it)) {
258 
259  Mat_Thermal temp_data;
260  CHKERR it->getAttributeDataStructure(temp_data);
261  setOfBlocks[it->getMeshsetId()].cOnductivity =
262  temp_data.data.Conductivity;
263  setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
264  CHKERR mField.get_moab().get_entities_by_type(
265  it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
267  setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
268 
269  Range skeleton;
270  CHKERR mField.get_moab().get_adjacencies(
271  setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
272  moab::Interface::UNION);
274  "MIX_SKELETON");
275  }
276 
277  // Define element to integrate natural boundary conditions, i.e. set values.
278  CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
280  fluxes_name);
282  fluxes_name);
284  fluxes_name);
286  values_name);
287 
288  // Define element to apply essential boundary conditions.
289  CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
291  fluxes_name);
293  fluxes_name);
295  fluxes_name);
297  values_name);
299  }
@ MF_ZERO
Definition: definitions.h:111
@ MAT_THERMALSET
block name is "MAT_THERMAL"
Definition: definitions.h:174
@ BLOCKSET
Definition: definitions.h:161
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
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 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 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 _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData

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

306  {
308  // build field
310  // get tetrahedrons which has been build previously and now in so called
311  // garbage bit level
312  Range done_tets;
313  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
314  BitRefLevel().set(0), BitRefLevel().set(), MBTET, done_tets);
315  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
316  BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTET,
317  done_tets);
318  // get tetrahedrons which belong to problem bit level
319  Range ref_tets;
320  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
321  ref_level, BitRefLevel().set(), MBTET, ref_tets);
322  ref_tets = subtract(ref_tets, done_tets);
323  CHKERR mField.build_finite_elements("MIX", &ref_tets, 2);
324  // get triangles which has been build previously and now in so called
325  // garbage bit level
326  Range done_faces;
327  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
328  BitRefLevel().set(0), BitRefLevel().set(), MBTRI, done_faces);
329  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
330  BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTRI,
331  done_faces);
332  // get triangles which belong to problem bit level
333  Range ref_faces;
334  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
335  ref_level, BitRefLevel().set(), MBTRI, ref_faces);
336  ref_faces = subtract(ref_faces, done_faces);
337  // build finite elements structures
338  CHKERR mField.build_finite_elements("MIX_BCFLUX", &ref_faces, 2);
339  CHKERR mField.build_finite_elements("MIX_BCVALUE", &ref_faces, 2);
340  CHKERR mField.build_finite_elements("MIX_SKELETON", &ref_faces, 2);
341  // Build adjacencies of degrees of freedom and elements
342  CHKERR mField.build_adjacencies(ref_level);
343  // Define problem
345  // set refinement level for problem
347  // Add element to problem
349  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_SKELETON");
350  // Boundary conditions
351  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCFLUX");
352  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCVALUE");
353  // build problem
354 
355  ProblemsManager *prb_mng_ptr;
356  CHKERR mField.getInterface(prb_mng_ptr);
357  CHKERR prb_mng_ptr->buildProblem("MIX", true);
358  // mesh partitioning
359  // partition
360  CHKERR prb_mng_ptr->partitionProblem("MIX");
361  CHKERR prb_mng_ptr->partitionFiniteElements("MIX");
362  // what are ghost nodes, see Petsc Manual
363  CHKERR prb_mng_ptr->partitionGhostDofs("MIX");
365  }
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:232
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ calculateResidual()

MoFEMErrorCode MixTransport::MixTransportElement::calculateResidual ( )

calculate residual

Definition at line 601 of file MixTransportElement.hpp.

601  {
603  CHKERR VecZeroEntries(F);
604  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
605  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
606  CHKERR VecAssemblyBegin(F);
607  CHKERR VecAssemblyEnd(F);
608  // calculate residuals
609  feVol.getOpPtrVector().clear();
610  feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
611  feVol.getOpPtrVector().push_back(
612  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
613  feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
614  feVol.getOpPtrVector().push_back(
615  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
616  feVol.getOpPtrVector().push_back(
617  new OpTauDotSigma_HdivHdiv(*this, "FLUXES", PETSC_NULL, F));
618  feVol.getOpPtrVector().push_back(
619  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", PETSC_NULL, F));
620  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
621  feTri.getOpPtrVector().clear();
622  feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
623  CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
624  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
625  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
626  CHKERR VecAssemblyBegin(F);
627  CHKERR VecAssemblyEnd(F);
628  // CHKERR VecAXPY(F,-1.,D0);
629  // CHKERR MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULL,F);
630  {
631  std::vector<int> ids;
632  ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
633  std::vector<double> vals(ids.size(), 0);
634  CHKERR VecSetValues(F, ids.size(), &*ids.begin(), &*vals.begin(),
635  INSERT_VALUES);
636  CHKERR VecAssemblyBegin(F);
637  CHKERR VecAssemblyEnd(F);
638  }
639  {
640  double nrm2_F;
641  CHKERR VecNorm(F, NORM_2, &nrm2_F);
642  PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
643  const double eps = 1e-8;
644  if (nrm2_F > eps) {
645  // SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"problem with
646  // residual");
647  }
648  }
650  }
static const double eps
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.

◆ createMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::createMatrices ( )

create matrices

Definition at line 462 of file MixTransportElement.hpp.

462  {
464  CHKERR mField.getInterface<MatrixManager>()
465  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("MIX", &Aij);
466  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
467  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
468  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
470  }
@ COL
Definition: definitions.h:136
@ ROW
Definition: definitions.h:136

◆ destroyMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::destroyMatrices ( )

destroy matrices

Definition at line 687 of file MixTransportElement.hpp.

687  {
689  CHKERR MatDestroy(&Aij);
690  ;
691  CHKERR VecDestroy(&D);
692  ;
693  CHKERR VecDestroy(&D0);
694  ;
695  CHKERR VecDestroy(&F);
696  ;
698  }

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

658  {
660  errorMap.clear();
661  sumErrorFlux = 0;
662  sumErrorDiv = 0;
663  sumErrorJump = 0;
664  feTri.getOpPtrVector().clear();
665  feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
666  CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
668  feVol.getOpPtrVector().clear();
669  feVol.getOpPtrVector().push_back(
670  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
671  feVol.getOpPtrVector().push_back(
672  new OpValuesGradientAtGaussPts(*this, "VALUES"));
673  feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
674  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
676  const Problem *problem_ptr;
677  CHKERR mField.get_problem("MIX", &problem_ptr);
678  PetscPrintf(mField.get_comm(),
679  "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
680  "jump^2 = %6.4e error tot^2 = %6.4e\n",
681  problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
684  }
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
map< double, EntityHandle > errorMap
virtual int get_comm_size() const =0
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 MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
double flux
impulse magnitude

◆ 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  }

◆ 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  }

◆ 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

◆ 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  }

◆ postProc()

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

Post process results.

Returns
error code

Definition at line 443 of file MixTransportElement.hpp.

443  {
446  CHKERR post_proc.generateReferenceElementMesh();
447  CHKERR post_proc.addFieldValuesPostProc("VALUES");
448  CHKERR post_proc.addFieldValuesGradientPostProc("VALUES");
449  CHKERR post_proc.addFieldValuesPostProc("FLUXES");
450  // CHKERR post_proc.addHdivFunctionsPostProc("FLUXES");
451  post_proc.getOpPtrVector().push_back(
452  new OpPostProc(post_proc.postProcMesh, post_proc.mapGaussPts));
453  CHKERR mField.loop_finite_elements("MIX", "MIX", post_proc);
454  CHKERR post_proc.writeFile(out_file.c_str());
456  }
Post processing.

◆ solveLinearProblem()

MoFEMErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )

solve problem

Returns
error code

Definition at line 476 of file MixTransportElement.hpp.

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

Member Data Documentation

◆ Aij

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

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

◆ D0

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

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

◆ F

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

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

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1348 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

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