v0.13.2
Loading...
Searching...
No Matches
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 29 of file MixTransportElement.hpp.

Constructor & Destructor Documentation

◆ MixTransportElement()

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

construction of this data structure

Definition at line 65 of file MixTransportElement.hpp.

66 : 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 ( )
inlinevirtual

destructor

Definition at line 71 of file MixTransportElement.hpp.

71{}

Member Function Documentation

◆ addFields()

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

Add fields to database.

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

Definition at line 189 of file MixTransportElement.hpp.

190 {
192 // Fields
195
196 // meshset consisting all entities in mesh
197 EntityHandle root_set = mField.get_moab().get_root_set();
198 // add entities to field
199 CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, fluxes);
200 CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, values);
201 CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
202 CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
203 CHKERR mField.set_field_order(root_set, MBTET, values, order);
205 }
@ 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 
)
inline

add finite elements

Definition at line 208 of file MixTransportElement.hpp.

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

Build problem.

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

Definition at line 291 of file MixTransportElement.hpp.

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

calculate residual

Definition at line 612 of file MixTransportElement.hpp.

612 {
614 CHKERR VecZeroEntries(F);
615 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
616 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
617 CHKERR VecAssemblyBegin(F);
618 CHKERR VecAssemblyEnd(F);
619 // calculate residuals
620 feVol.getOpPtrVector().clear();
621 CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
622 feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
623 feVol.getOpPtrVector().push_back(
624 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
625 feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
626 feVol.getOpPtrVector().push_back(
627 new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
628 feVol.getOpPtrVector().push_back(
629 new OpTauDotSigma_HdivHdiv(*this, "FLUXES", PETSC_NULL, F));
630 feVol.getOpPtrVector().push_back(
631 new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", PETSC_NULL, F));
633 feTri.getOpPtrVector().clear();
634 feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
635 CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
636 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
637 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
638 CHKERR VecAssemblyBegin(F);
639 CHKERR VecAssemblyEnd(F);
640 // CHKERR VecAXPY(F,-1.,D0);
641 // CHKERR MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULL,F);
642 {
643 std::vector<int> ids;
644 ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
645 std::vector<double> vals(ids.size(), 0);
646 CHKERR VecSetValues(F, ids.size(), &*ids.begin(), &*vals.begin(),
647 INSERT_VALUES);
648 CHKERR VecAssemblyBegin(F);
649 CHKERR VecAssemblyEnd(F);
650 }
651 {
652 double nrm2_F;
653 CHKERR VecNorm(F, NORM_2, &nrm2_F);
654 PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
655 const double eps = 1e-8;
656 if (nrm2_F > eps) {
657 // SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"problem with
658 // residual");
659 }
660 }
662 }
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 ( )
inline

create matrices

Definition at line 474 of file MixTransportElement.hpp.

474 {
476 CHKERR mField.getInterface<MatrixManager>()
477 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("MIX", &Aij);
478 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
479 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
480 CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
482 }
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123

◆ destroyMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::destroyMatrices ( )
inline

destroy matrices

Definition at line 700 of file MixTransportElement.hpp.

700 {
702 CHKERR MatDestroy(&Aij);
703 ;
704 CHKERR VecDestroy(&D);
705 ;
706 CHKERR VecDestroy(&D0);
707 ;
708 CHKERR VecDestroy(&F);
709 ;
711 }

◆ evaluateError()

MoFEMErrorCode MixTransport::MixTransportElement::evaluateError ( )
inline

Calculate error on elements.

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

Definition at line 670 of file MixTransportElement.hpp.

670 {
672 errorMap.clear();
673 sumErrorFlux = 0;
674 sumErrorDiv = 0;
675 sumErrorJump = 0;
676 feTri.getOpPtrVector().clear();
677 feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
678 CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
680 feVol.getOpPtrVector().clear();
681 CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
682 feVol.getOpPtrVector().push_back(
683 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
684 feVol.getOpPtrVector().push_back(
685 new OpValuesGradientAtGaussPts(*this, "VALUES"));
686 feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
687 CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
689 const Problem *problem_ptr;
690 CHKERR mField.get_problem("MIX", &problem_ptr);
691 PetscPrintf(mField.get_comm(),
692 "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
693 "jump^2 = %6.4e error tot^2 = %6.4e\n",
694 problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
697 }
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 
)
inlinevirtual

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

165 {
167 flux = 0;
169 }
#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 
)
inlinevirtual

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

148 {
150 value = 0;
152 }

◆ getDirichletBCIndices()

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

get dof indices where essential boundary conditions are applied

Parameters
isindices
Returns
error code

Definition at line 87 of file MixTransportElement.hpp.

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

◆ getResistivity()

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

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

128 {
130 inv_k.clear();
131 for (int dd = 0; dd < 3; dd++) {
132 inv_k(dd, dd) = 1;
133 }
135 }
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 
)
inlinevirtual

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

111 {
113 flux = 0;
115 }

◆ postProc()

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

Post process results.

Returns
error code

Definition at line 428 of file MixTransportElement.hpp.

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

◆ solveLinearProblem()

MoFEMErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )
inline

solve problem

Returns
error code

Definition at line 488 of file MixTransportElement.hpp.

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

Member Data Documentation

◆ Aij

Mat MixTransport::MixTransportElement::Aij

Definition at line 471 of file MixTransportElement.hpp.

◆ bcIndices

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

Definition at line 80 of file MixTransportElement.hpp.

◆ D

Vec MixTransport::MixTransportElement::D

Definition at line 470 of file MixTransportElement.hpp.

◆ D0

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

Definition at line 470 of file MixTransportElement.hpp.

◆ divergenceAtGaussPts

VectorDouble MixTransport::MixTransportElement::divergenceAtGaussPts

divergence at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 77 of file MixTransportElement.hpp.

◆ errorMap

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

Definition at line 1360 of file MixTransportElement.hpp.

◆ F

Vec MixTransport::MixTransportElement::F

Definition at line 470 of file MixTransportElement.hpp.

◆ feTri

MyTriFE MixTransport::MixTransportElement::feTri

Instance of surface element.

Definition at line 60 of file MixTransportElement.hpp.

◆ feVol

MyVolumeFE MixTransport::MixTransportElement::feVol

Instance of volume element.

Definition at line 46 of file MixTransportElement.hpp.

◆ fluxesAtGaussPts

MatrixDouble MixTransport::MixTransportElement::fluxesAtGaussPts

fluxes at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 78 of file MixTransportElement.hpp.

◆ mField

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

Definition at line 31 of file MixTransportElement.hpp.

◆ setOfBlocks

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

maps block set id with appropriate BlockData

Definition at line 180 of file MixTransportElement.hpp.

◆ sumErrorDiv

double MixTransport::MixTransportElement::sumErrorDiv

Definition at line 1362 of file MixTransportElement.hpp.

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1361 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

Definition at line 1363 of file MixTransportElement.hpp.

◆ valuesAtGaussPts

VectorDouble MixTransport::MixTransportElement::valuesAtGaussPts

values at integration points on element

Examples
UnsaturatedFlow.hpp.

Definition at line 73 of file MixTransportElement.hpp.

◆ valuesGradientAtGaussPts

MatrixDouble MixTransport::MixTransportElement::valuesGradientAtGaussPts

gradients at integration points on element

Definition at line 75 of file MixTransportElement.hpp.


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