v0.13.1
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, EntityHandle > errorMap
 
double sumErrorFlux
 
double sumErrorDiv
 
double sumErrorJump
 

Detailed Description

Mix transport problem.

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

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

General problem have form,

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

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

68 : 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 73 of file MixTransportElement.hpp.

73{}

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

192 {
194 // Fields
197
198 // meshset consisting all entities in mesh
199 EntityHandle root_set = mField.get_moab().get_root_set();
200 // add entities to field
201 CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, fluxes);
202 CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, values);
203 CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
204 CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
205 CHKERR mField.set_field_order(root_set, MBTET, values, order);
207 }
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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 moab::Interface & get_moab()=0
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.

◆ addFiniteElements()

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

add finite elements

Definition at line 210 of file MixTransportElement.hpp.

211 {
213
214 // Set up volume element operators. Operators are used to calculate
215 // components of stiffness matrix & right hand side, in essence are used to
216 // do volume integrals over tetrahedral in this case.
217
218 // Define element "MIX". Note that this element will work with fluxes_name
219 // and values_name. This reflect bilinear form for the problem
227
228 // Define finite element to integrate over skeleton, we need that to
229 // evaluate error
230 CHKERR mField.add_finite_element("MIX_SKELETON", MF_ZERO);
232 fluxes_name);
234 fluxes_name);
236 fluxes_name);
237
238 // Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those
239 // BLOCKSETS and add them to "MIX" finite element. In addition get data form
240 // that meshset and set conductivity which is used to calculate fluxes from
241 // gradients of concentration or gradient of temperature, depending how you
242 // interpret variables.
245
246 Mat_Thermal temp_data;
247 CHKERR it->getAttributeDataStructure(temp_data);
248 setOfBlocks[it->getMeshsetId()].cOnductivity =
249 temp_data.data.Conductivity;
250 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
251 CHKERR mField.get_moab().get_entities_by_type(
252 it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
254 setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
255
256 Range skeleton;
257 CHKERR mField.get_moab().get_adjacencies(
258 setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
259 moab::Interface::UNION);
261 "MIX_SKELETON");
262 }
263
264 // Define element to integrate natural boundary conditions, i.e. set values.
265 CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
267 fluxes_name);
269 fluxes_name);
271 fluxes_name);
273 values_name);
274
275 // Define element to apply essential boundary conditions.
278 fluxes_name);
280 fluxes_name);
282 fluxes_name);
284 values_name);
286 }
@ MF_ZERO
Definition: definitions.h:98
@ MAT_THERMALSET
block name is "MAT_THERMAL"
Definition: definitions.h:161
@ BLOCKSET
Definition: definitions.h:148
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 293 of file MixTransportElement.hpp.

293 {
295 // build field
297 // get tetrahedrons which has been build previously and now in so called
298 // garbage bit level
299 Range done_tets;
300 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
301 BitRefLevel().set(0), BitRefLevel().set(), MBTET, done_tets);
302 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
303 BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTET,
304 done_tets);
305 // get tetrahedrons which belong to problem bit level
306 Range ref_tets;
307 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
308 ref_level, BitRefLevel().set(), MBTET, ref_tets);
309 ref_tets = subtract(ref_tets, done_tets);
310 CHKERR mField.build_finite_elements("MIX", &ref_tets, 2);
311 // get triangles which has been build previously and now in so called
312 // garbage bit level
313 Range done_faces;
314 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
315 BitRefLevel().set(0), BitRefLevel().set(), MBTRI, done_faces);
316 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
317 BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTRI,
318 done_faces);
319 // get triangles which belong to problem bit level
320 Range ref_faces;
321 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
322 ref_level, BitRefLevel().set(), MBTRI, ref_faces);
323 ref_faces = subtract(ref_faces, done_faces);
324 // build finite elements structures
325 CHKERR mField.build_finite_elements("MIX_BCFLUX", &ref_faces, 2);
326 CHKERR mField.build_finite_elements("MIX_BCVALUE", &ref_faces, 2);
327 CHKERR mField.build_finite_elements("MIX_SKELETON", &ref_faces, 2);
328 // Build adjacencies of degrees of freedom and elements
330 // Define problem
332 // set refinement level for problem
334 // Add element to problem
336 CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_SKELETON");
337 // Boundary conditions
339 CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCVALUE");
340 // build problem
341
342 ProblemsManager *prb_mng_ptr;
343 CHKERR mField.getInterface(prb_mng_ptr);
344 CHKERR prb_mng_ptr->buildProblem("MIX", true);
345 // mesh partitioning
346 // partition
347 CHKERR prb_mng_ptr->partitionProblem("MIX");
348 CHKERR prb_mng_ptr->partitionFiniteElements("MIX");
349 // what are ghost nodes, see Petsc Manual
350 CHKERR prb_mng_ptr->partitionGhostDofs("MIX");
352 }
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
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:40
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 624 of file MixTransportElement.hpp.

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

485 {
487 CHKERR mField.getInterface<MatrixManager>()
488 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("MIX", &Aij);
489 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
490 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
491 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
493 }
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123

◆ destroyMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::destroyMatrices ( )

destroy matrices

Definition at line 710 of file MixTransportElement.hpp.

710 {
712 CHKERR MatDestroy(&Aij);
713 ;
714 CHKERR VecDestroy(&D);
715 ;
716 CHKERR VecDestroy(&D0);
717 ;
718 CHKERR VecDestroy(&F);
719 ;
721 }

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

681 {
683 errorMap.clear();
684 sumErrorFlux = 0;
685 sumErrorDiv = 0;
686 sumErrorJump = 0;
687 feTri.getOpPtrVector().clear();
688 feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
689 CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
691 feVol.getOpPtrVector().clear();
692 feVol.getOpPtrVector().push_back(
693 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
694 feVol.getOpPtrVector().push_back(
695 new OpValuesGradientAtGaussPts(*this, "VALUES"));
696 feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
697 CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
699 const Problem *problem_ptr;
700 CHKERR mField.get_problem("MIX", &problem_ptr);
701 PetscPrintf(mField.get_comm(),
702 "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
703 "jump^2 = %6.4e error tot^2 = %6.4e\n",
704 problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
707 }
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 MyTransport, MyTransport, and MixTransport::UnsaturatedFlowElement.

Definition at line 165 of file MixTransportElement.hpp.

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

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

150 {
152 value = 0;
154 }

◆ getDirichletBCIndices()

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

get dof indices where essential boundary conditions are applied

Parameters
isindices
Returns
error code

Definition at line 89 of file MixTransportElement.hpp.

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

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

130 {
132 inv_k.clear();
133 for (int dd = 0; dd < 3; dd++) {
134 inv_k(dd, dd) = 1;
135 }
137 }
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 111 of file MixTransportElement.hpp.

113 {
115 flux = 0;
117 }

◆ postProc()

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

Post process results.

Returns
error code

Definition at line 430 of file MixTransportElement.hpp.

430 {
432
433 PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
434 mField);
435
436 auto jac_ptr = boost::make_shared<MatrixDouble>();
437 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
438 auto det_ptr = boost::make_shared<VectorDouble>();
439 post_proc.getOpPtrVector().push_back(new OpCalculateHOJac<3>(jac_ptr));
440 post_proc.getOpPtrVector().push_back(
441 new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
442 post_proc.getOpPtrVector().push_back(
443 new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
444 post_proc.getOpPtrVector().push_back(
445 new OpSetHOInvJacToScalarBases<3>(L2, inv_jac_ptr));
446
447 auto values_ptr = boost::make_shared<VectorDouble>();
448 auto grad_ptr = boost::make_shared<MatrixDouble>();
449 auto flux_ptr = boost::make_shared<MatrixDouble>();
450
451 post_proc.getOpPtrVector().push_back(
452 new OpCalculateScalarFieldValues("VALUES", values_ptr));
453 post_proc.getOpPtrVector().push_back(
454 new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
455 post_proc.getOpPtrVector().push_back(
456 new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
457
458 using OpPPMap = OpPostProcMapInMoab<3, 3>;
459
460 post_proc.getOpPtrVector().push_back(
461
462 new OpPPMap(post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
463
464 {{"VALUES", values_ptr}},
465
466 {{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
467
468 {}, {}
469
470 ));
471
472 post_proc.getOpPtrVector().push_back(new OpPostProc(
473 post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
474
475 CHKERR mField.loop_finite_elements("MIX", "MIX", post_proc);
476
477 CHKERR post_proc.writeFile(out_file.c_str());
479 }
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Post post-proc data at points from hash maps.

◆ solveLinearProblem()

MoFEMErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )

solve problem

Returns
error code

Definition at line 499 of file MixTransportElement.hpp.

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

Member Data Documentation

◆ Aij

Mat MixTransport::MixTransportElement::Aij

Definition at line 482 of file MixTransportElement.hpp.

◆ bcIndices

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

Definition at line 82 of file MixTransportElement.hpp.

◆ D

Vec MixTransport::MixTransportElement::D

Definition at line 481 of file MixTransportElement.hpp.

◆ D0

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

Definition at line 481 of file MixTransportElement.hpp.

◆ divergenceAtGaussPts

VectorDouble MixTransport::MixTransportElement::divergenceAtGaussPts

divergence at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 79 of file MixTransportElement.hpp.

◆ errorMap

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

Definition at line 1370 of file MixTransportElement.hpp.

◆ F

Vec MixTransport::MixTransportElement::F

Definition at line 481 of file MixTransportElement.hpp.

◆ feTri

MyTriFE MixTransport::MixTransportElement::feTri

Instance of surface element.

Definition at line 62 of file MixTransportElement.hpp.

◆ feVol

MyVolumeFE MixTransport::MixTransportElement::feVol

Instance of volume element.

Definition at line 48 of file MixTransportElement.hpp.

◆ fluxesAtGaussPts

MatrixDouble MixTransport::MixTransportElement::fluxesAtGaussPts

fluxes at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 80 of file MixTransportElement.hpp.

◆ mField

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

Definition at line 33 of file MixTransportElement.hpp.

◆ setOfBlocks

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

maps block set id with appropriate BlockData

Definition at line 182 of file MixTransportElement.hpp.

◆ sumErrorDiv

double MixTransport::MixTransportElement::sumErrorDiv

Definition at line 1372 of file MixTransportElement.hpp.

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1371 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

Definition at line 1373 of file MixTransportElement.hpp.

◆ valuesAtGaussPts

VectorDouble MixTransport::MixTransportElement::valuesAtGaussPts

values at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 75 of file MixTransportElement.hpp.

◆ valuesGradientAtGaussPts

MatrixDouble MixTransport::MixTransportElement::valuesGradientAtGaussPts

gradients at integration points on element

Definition at line 77 of file MixTransportElement.hpp.


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