Mix transport problemNote to solve this system you need to use direct solver or proper preconditioner for saddle problem. More...

#include <users_modules/basic_finite_elements/mix_transport/src/MixTransportElement.hpp>

Inheritance diagram for MixTransport::MixTransportElement:

[legend]

Collaboration diagram for MixTransport::MixTransportElement:

[legend]

Classes
struct	BlockData
	data for evaluation of het conductivity and heat capacity elements mofem_thermal_elem 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...

PetscErrorCode	getDirichletBCIndices (IS *is)
	get dof indices where essential boundary conditions are applied More...

virtual PetscErrorCode	getSource (const EntityHandle ent, const double x, const double y, const double z, double &flux)
	set source term More...

virtual PetscErrorCode	getResistivity (const EntityHandle ent, const double x, const double y, const double z, MatrixDouble3by3 &inv_k)
	natural (Dirihlet) boundary conditions (set values) More...

virtual PetscErrorCode	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 PetscErrorCode	getBcOnFluxes (const EntityHandle ent, const double x, const double y, const double z, double &flux)
	essential (Neumann) boundary condition (set fluxes) More...

PetscErrorCode	addFields (const std::string &values, const std::string &fluxes, const int order)
	Add fields to database. More...

PetscErrorCode	addFiniteElements (const std::string &fluxes_name, const std::string &values_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
	add finite elements More...

PetscErrorCode	buildProblem (BitRefLevel &ref_level)
	Build problem. More...

PetscErrorCode	postProc (const string out_file)
	Post process results. More...

PetscErrorCode	createMatrices ()
	create matrices More...

PetscErrorCode	solveLinearProblem ()
	solve problem More...

PetscErrorCode	calculateResidual ()
	calculate residual More...

PetscErrorCode	evaluateError ()
	Calculate error on elements. More...

PetscErrorCode	destroyMatrices ()
	destroy matrices More...

Public Attributes
MoFEM::Interface &	mField

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, BlockData >	setOfBlocks
	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 [5] [6] [16] https://www.researchgate.net/profile/Richard_Falk/publication/226454406_Differential_Complexes_and_Stability_of_Finite_Element_Methods_I._The_de_Rham_Complex/links/02e7e5214f0426ff77000000.pdf

General problem have form,

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

Definition at line 44 of file MixTransportElement.hpp.

Constructor & Destructor Documentation

◆ MixTransportElement()

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

construction of this data structure

Definition at line 78 of file MixTransportElement.hpp.

                                               :
   mField(m_field),
   feVol(m_field),
   feTri(m_field) {};

◆ ~MixTransportElement()

virtual MixTransport::MixTransportElement::~MixTransportElement ( )

virtual

destructor

Definition at line 86 of file MixTransportElement.hpp.

86 {}

Member Function Documentation

◆ addFields()

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

Add fields to database.

Parameters

values	name of the fields
fluxes	name of filed for fluxes
order	order of approximation

Returns: error code

Definition at line 208 of file MixTransportElement.hpp.

                                                                                             {
 
     PetscFunctionBegin;
     //Fields
     ierr = mField.add_field(fluxes,HDIV,DEMKOWICZ_JACOBI_BASE,1); CHKERRQ(ierr);
     ierr = mField.add_field(values,L2,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
 
     //meshset consisting all entities in mesh
     EntityHandle root_set = mField.get_moab().get_root_set();
     //add entities to field
     ierr = mField.add_ents_to_field_by_type(root_set,MBTET,fluxes); CHKERRQ(ierr);
     ierr = mField.add_ents_to_field_by_type(root_set,MBTET,values); CHKERRQ(ierr);
     ierr = mField.set_field_order(root_set,MBTET,fluxes,order+1); CHKERRQ(ierr);
     ierr = mField.set_field_order(root_set,MBTRI,fluxes,order+1); CHKERRQ(ierr);
     ierr = mField.set_field_order(root_set,MBTET,values,order); CHKERRQ(ierr);
 
 
     PetscFunctionReturn(0);
   }

◆ addFiniteElements()

PetscErrorCode MixTransport::MixTransportElement::addFiniteElements	(	const std::string &	fluxes_name,
		const std::string &	values_name,
		const std::string	mesh_nodals_positions = `"MESH_NODE_POSITIONS"`
	)

add finite elements

Definition at line 229 of file MixTransportElement.hpp.

     {
     PetscFunctionBegin;
 
 
 
 
     // Set up volume element operators. Operators are used to calculate components
     // of stiffness matrix & right hand side, in essence are used to do volume integrals over
     // tetrahedral in this case.
 
     // Define element "MIX". Note that this element will work with fluxes_name and
     // values_name. This reflect bilinear form for the problem
     ierr = mField.add_finite_element("MIX",MF_ZERO); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_row("MIX",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_col("MIX",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_row("MIX",values_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_col("MIX",values_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_data("MIX",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_data("MIX",values_name); CHKERRQ(ierr);
 
     // Define finite element to integrate over skeleton, we need that to evaluate error
     ierr = mField.add_finite_element("MIX_SKELETON",MF_ZERO); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_row("MIX_SKELETON",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_col("MIX_SKELETON",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_data("MIX_SKELETON",fluxes_name); CHKERRQ(ierr);
 
     // In some cases you like to use HO geometry to describe shape of the bode, curved edges and faces, for
     // example body is a sphere. HO geometry is approximated by a field,  which can be hierarchical, so shape of
     // the edges could be given by polynomial of arbitrary order.
     //
     // Check if field "mesh_nodals_positions" is defined, and if it is add that field to data of finite
     // element. MoFEM will use that that to calculate Jacobian as result that geometry in nonlinear.
     if(mField.check_field(mesh_nodals_positions)) {
       ierr = mField.modify_finite_element_add_field_data("MIX",mesh_nodals_positions); CHKERRQ(ierr);
       ierr = mField.modify_finite_element_add_field_data("MIX_SKELETON",mesh_nodals_positions); CHKERRQ(ierr);
     }
     // Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those BLOCKSETS
     // and add them to "MIX" finite element. In addition get data form that meshset
     // and set cOnductivity which is used to calculate fluxes from gradients of concentration
     // or gradient of temperature, depending how you interpret variables.
     for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField,BLOCKSET|MAT_THERMALSET,it)) {
 
       // cerr << *it << endl;
 
       Mat_Thermal temp_data;
       ierr = it->getAttributeDataStructure(temp_data); CHKERRQ(ierr);
       setOfBlocks[it->getMeshsetId()].cOnductivity = temp_data.data.Conductivity;
       setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
       rval = mField.get_moab().get_entities_by_type(
         it->meshset,MBTET,setOfBlocks[it->getMeshsetId()].tEts,true
       ); CHKERRQ_MOAB(rval);
       ierr = mField.add_ents_to_finite_element_by_type(
         setOfBlocks[it->getMeshsetId()].tEts,MBTET,"MIX"
       ); CHKERRQ(ierr);
 
       Range skeleton;
       rval = mField.get_moab().get_adjacencies(
         setOfBlocks[it->getMeshsetId()].tEts,2,false,skeleton,moab::Interface::UNION
       );
       ierr = mField.add_ents_to_finite_element_by_type(
         skeleton,MBTRI,"MIX_SKELETON"
       ); CHKERRQ(ierr);
 
     }
 
     // Define element to integrate natural boundary conditions, i.e. set values.
     ierr = mField.add_finite_element("MIX_BCVALUE",MF_ZERO); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_row("MIX_BCVALUE",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_col("MIX_BCVALUE",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_data("MIX_BCVALUE",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_data("MIX_BCVALUE",values_name); CHKERRQ(ierr);
     if(mField.check_field(mesh_nodals_positions)) {
       ierr = mField.modify_finite_element_add_field_data("MIX_BCVALUE",mesh_nodals_positions); CHKERRQ(ierr);
     }
 
     // Define element to apply essential boundary conditions.
     ierr = mField.add_finite_element("MIX_BCFLUX",MF_ZERO); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_row("MIX_BCFLUX",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_col("MIX_BCFLUX",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_data("MIX_BCFLUX",fluxes_name); CHKERRQ(ierr);
     ierr = mField.modify_finite_element_add_field_data("MIX_BCFLUX",values_name); CHKERRQ(ierr);
     if(mField.check_field(mesh_nodals_positions)) {
       ierr = mField.modify_finite_element_add_field_data("MIX_BCFLUX",mesh_nodals_positions); CHKERRQ(ierr);
     }
 
     PetscFunctionReturn(0);
   }

◆ buildProblem()

PetscErrorCode MixTransport::MixTransportElement::buildProblem ( BitRefLevel & ref_level )

Build problem.

Parameters

ref_level mesh refinement on which mesh problem you like to built.

Returns: error code

Definition at line 326 of file MixTransportElement.hpp.

                                                       {
 
     PetscFunctionBegin;
     //build field
     ierr = mField.build_fields(); CHKERRQ(ierr);
     // get tetrahedrons which has been build previously and now in so called garbage bit level
     Range done_tets;
     ierr = mField.get_entities_by_type_and_ref_level(
       BitRefLevel().set(0),BitRefLevel().set(),MBTET,done_tets
     ); CHKERRQ(ierr);
     ierr = mField.get_entities_by_type_and_ref_level(
       BitRefLevel().set(BITREFLEVEL_SIZE-1),BitRefLevel().set(),MBTET,done_tets
     ); CHKERRQ(ierr);
     // get tetrahedrons which belong to problem bit level
     Range ref_tets;
     ierr = mField.get_entities_by_type_and_ref_level(
       ref_level,BitRefLevel().set(),MBTET,ref_tets
     ); CHKERRQ(ierr);
     ref_tets = subtract(ref_tets,done_tets);
     ierr = mField.build_finite_elements("MIX",&ref_tets,2); CHKERRQ(ierr);
     // get triangles which has been build previously and now in so called garbage bit level
     Range done_faces;
     ierr = mField.get_entities_by_type_and_ref_level(
       BitRefLevel().set(0),BitRefLevel().set(),MBTRI,done_faces
     ); CHKERRQ(ierr);
     ierr = mField.get_entities_by_type_and_ref_level(
       BitRefLevel().set(BITREFLEVEL_SIZE-1),BitRefLevel().set(),MBTRI,done_faces
     ); CHKERRQ(ierr);
     // get triangles which belong to problem bit level
     Range ref_faces;
     ierr = mField.get_entities_by_type_and_ref_level(
       ref_level,BitRefLevel().set(),MBTRI,ref_faces
     ); CHKERRQ(ierr);
     ref_faces = subtract(ref_faces,done_faces);
     //build finite elements structures
     ierr = mField.build_finite_elements("MIX_BCFLUX",&ref_faces,2); CHKERRQ(ierr);
     ierr = mField.build_finite_elements("MIX_BCVALUE",&ref_faces,2); CHKERRQ(ierr);
     ierr = mField.build_finite_elements("MIX_SKELETON",&ref_faces,2); CHKERRQ(ierr);
     //Build adjacencies of degrees of freedom and elements
     ierr = mField.build_adjacencies(ref_level); CHKERRQ(ierr);
     //Define problem
     ierr = mField.add_problem("MIX",MF_ZERO); CHKERRQ(ierr);
     //set refinement level for problem
     ierr = mField.modify_problem_ref_level_set_bit("MIX",ref_level); CHKERRQ(ierr);
     // Add element to problem
     ierr = mField.modify_problem_add_finite_element("MIX","MIX"); CHKERRQ(ierr);
     ierr = mField.modify_problem_add_finite_element("MIX","MIX_SKELETON"); CHKERRQ(ierr);
     // Boundary conditions
     ierr = mField.modify_problem_add_finite_element("MIX","MIX_BCFLUX"); CHKERRQ(ierr);
     ierr = mField.modify_problem_add_finite_element("MIX","MIX_BCVALUE"); CHKERRQ(ierr);
     //build problem
 
     ProblemsManager *prb_mng_ptr;
     ierr = mField.query_interface(prb_mng_ptr); CHKERRQ(ierr);
     ierr = prb_mng_ptr->buildProblem("MIX",true); CHKERRQ(ierr);
     //mesh partitioning
     //partition
     ierr = prb_mng_ptr->partitionProblem("MIX"); CHKERRQ(ierr);
     ierr = prb_mng_ptr->partitionFiniteElements("MIX"); CHKERRQ(ierr);
     //what are ghost nodes, see Petsc Manual
     ierr = prb_mng_ptr->partitionGhostDofs("MIX"); CHKERRQ(ierr);
     PetscFunctionReturn(0);
   }

◆ calculateResidual()

PetscErrorCode MixTransport::MixTransportElement::calculateResidual ( )

calculate residual

Definition at line 613 of file MixTransportElement.hpp.

                                      {
 
     PetscFunctionBegin;
     ierr = VecZeroEntries(F); CHKERRQ(ierr);
     ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
     ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
     ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
     ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
     //calculate residuals
     feVol.getOpPtrVector().clear();
     feVol.getOpPtrVector().push_back(new OpL2Source(*this,"VALUES",F));
     feVol.getOpPtrVector().push_back(new OpFluxDivergenceAtGaussPts(*this,"FLUXES"));
     feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this,"VALUES"));
     feVol.getOpPtrVector().push_back(new OpDivTauU_HdivL2(*this,"FLUXES","VALUES",F));
     feVol.getOpPtrVector().push_back(new OpTauDotSigma_HdivHdiv(*this,"FLUXES",PETSC_NULL,F));
     feVol.getOpPtrVector().push_back(new OpVDivSigma_L2Hdiv(*this,"VALUES","FLUXES",PETSC_NULL,F));
     ierr = mField.loop_finite_elements("MIX","MIX",feVol); CHKERRQ(ierr);
     feTri.getOpPtrVector().clear();
     feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this,"FLUXES",F));
     ierr = mField.loop_finite_elements("MIX","MIX_BCVALUE",feTri); CHKERRQ(ierr);
     ierr = VecGhostUpdateBegin(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
     ierr = VecGhostUpdateEnd(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
     ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
     ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
     // ierr = VecAXPY(F,-1.,D0); CHKERRQ(ierr);
     // ierr = MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULL,F); CHKERRQ(ierr);
     {
       std::vector<int> ids;
       ids.insert(ids.begin(),bcIndices.begin(),bcIndices.end());
       std::vector<double> vals(ids.size(),0);
       ierr = VecSetValues(F,ids.size(),&*ids.begin(),&*vals.begin(),INSERT_VALUES); CHKERRQ(ierr);
       ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
       ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
     }
     {
       double nrm2_F;
       ierr = VecNorm(F,NORM_2,&nrm2_F); CHKERRQ(ierr);
       PetscPrintf(PETSC_COMM_WORLD,"nrm2_F = %6.4e\n",nrm2_F);
       const double eps = 1e-8;
       if(nrm2_F > eps) {
         //SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"problem with residual");
       }
     }
     PetscFunctionReturn(0);
   }

◆ createMatrices()

PetscErrorCode MixTransport::MixTransportElement::createMatrices ( )

create matrices

Definition at line 484 of file MixTransportElement.hpp.

                                   {
     PetscFunctionBegin;
     ierr = mField.MatCreateMPIAIJWithArrays("MIX",&Aij); CHKERRQ(ierr);
     ierr = mField.query_interface<VecManager>()->vecCreateGhost("MIX",COL,&D); CHKERRQ(ierr);
     ierr = mField.query_interface<VecManager>()->vecCreateGhost("MIX",COL,&D0); CHKERRQ(ierr);
     ierr = mField.query_interface<VecManager>()->vecCreateGhost("MIX",ROW,&F); CHKERRQ(ierr);
     PetscFunctionReturn(0);
   }

◆ destroyMatrices()

PetscErrorCode MixTransport::MixTransportElement::destroyMatrices ( )

destroy matrices

Definition at line 692 of file MixTransportElement.hpp.

                                    {
     PetscFunctionBegin;
     ierr = MatDestroy(&Aij); CHKERRQ(ierr);
     ierr = VecDestroy(&D); CHKERRQ(ierr);
     ierr = VecDestroy(&D0); CHKERRQ(ierr);
     ierr = VecDestroy(&F); CHKERRQ(ierr);
     PetscFunctionReturn(0);
   }

◆ evaluateError()

PetscErrorCode MixTransport::MixTransportElement::evaluateError ( )

Calculate error on elements.

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

Definition at line 665 of file MixTransportElement.hpp.

                                  {
 
     PetscFunctionBegin;
     errorMap.clear();
     sumErrorFlux = 0;
     sumErrorDiv = 0;
     sumErrorJump = 0;
     feTri.getOpPtrVector().clear();
     feTri.getOpPtrVector().push_back(new OpSkeleton(*this,"FLUXES"));
     ierr = mField.loop_finite_elements("MIX","MIX_SKELETON",feTri,0,mField.get_comm_size()); CHKERRQ(ierr);
     feVol.getOpPtrVector().clear();
     feVol.getOpPtrVector().push_back(new OpFluxDivergenceAtGaussPts(*this,"FLUXES"));
     feVol.getOpPtrVector().push_back(new OpValuesGradientAtGaussPts(*this,"VALUES"));
     feVol.getOpPtrVector().push_back(new OpError(*this,"VALUES"));
     ierr = mField.loop_finite_elements("MIX","MIX",feVol,0,mField.get_comm_size()); CHKERRQ(ierr);
     const Problem *problem_ptr;
     ierr = mField.get_problem("MIX",&problem_ptr); CHKERRQ(ierr);
     PetscPrintf(
       mField.get_comm(),
       "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error jump^2 = %6.4e error tot^2 = %6.4e\n",
       problem_ptr->getNbDofsRow(),
       sumErrorFlux,sumErrorDiv,sumErrorJump,sumErrorFlux+sumErrorDiv+sumErrorJump
     );
     PetscFunctionReturn(0);
   }

◆ getBcOnFluxes()

virtual PetscErrorCode 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

ent	handle to finite element entity
x	coord
y	coord
z	coord
flux	reference to flux which is set by function

Returns: [description]

Reimplemented in MixTransport::UnsaturatedFlowElement, MyTransport, and MyTransport.

Definition at line 182 of file MixTransportElement.hpp.

                   {
     PetscFunctionBegin;
     flux = 0;
     PetscFunctionReturn(0);
   }

◆ getBcOnValues()

virtual PetscErrorCode 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

ent	entity on which bc is evaluated
x	coordinate
y	coordinate
z	coordinate
value	vale

Returns: error code

Reimplemented in MixTransport::UnsaturatedFlowElement.

Definition at line 163 of file MixTransportElement.hpp.

                    {
     PetscFunctionBegin;
     value = 0;
     PetscFunctionReturn(0);
   }

◆ getDirichletBCIndices()

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

get dof indices where essential boundary conditions are applied

Parameters

is indices

Returns: error code

Definition at line 100 of file MixTransportElement.hpp.

                                                {
     PetscFunctionBegin;
     std::vector<int> ids;
     ids.insert(ids.begin(),bcIndices.begin(),bcIndices.end());
     IS is_local;
     ierr = ISCreateGeneral(
       mField.get_comm(),ids.size(),ids.empty()?PETSC_NULL:&ids[0],PETSC_COPY_VALUES,&is_local
     ); CHKERRQ(ierr);
     ierr = ISAllGather(is_local,is); CHKERRQ(ierr);
     ierr = ISDestroy(&is_local); CHKERRQ(ierr);
     PetscFunctionReturn(0);
   }

◆ getResistivity()

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

virtual

natural (Dirihlet) boundary conditions (set values)

Parameters

ent	handle to finite element entity
x	coord
y	coord
z	coord
value	reference to value set by function

Returns: error code

Definition at line 141 of file MixTransportElement.hpp.

     {
     PetscFunctionBegin;
     inv_k.clear();
     for(int dd = 0;dd<3;dd++) {
       inv_k(dd,dd) = 1;
     }
     PetscFunctionReturn(0);
   }

◆ getSource()

virtual PetscErrorCode MixTransport::MixTransportElement::getSource	(	const EntityHandle	ent,
		const double	x,
		const double	y,
		const double	z,
		double &	flux
	)

virtual

set source term

Parameters

ent	handle to entity on which function is evaluated
x	coord
y	coord
z	coord
flux	reference to source term set by function

Returns: error code

Reimplemented in MyTransport, and MyTransport.

Definition at line 122 of file MixTransportElement.hpp.

     {
     PetscFunctionBegin;
     flux  = 0;
     PetscFunctionReturn(0);
   }

◆ postProc()

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

Post process results.

Returns: error code

Definition at line 465 of file MixTransportElement.hpp.

                                                  {
 
     PetscFunctionBegin;
     PostProcVolumeOnRefinedMesh post_proc(mField);
     ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
     ierr = post_proc.addFieldValuesPostProc("VALUES"); CHKERRQ(ierr);
     ierr = post_proc.addFieldValuesGradientPostProc("VALUES"); CHKERRQ(ierr);
     ierr = post_proc.addFieldValuesPostProc("FLUXES"); CHKERRQ(ierr);
     // ierr = post_proc.addHdivFunctionsPostProc("FLUXES"); CHKERRQ(ierr);
     post_proc.getOpPtrVector().push_back(new OpPostProc(post_proc.postProcMesh,post_proc.mapGaussPts));
     ierr = mField.loop_finite_elements("MIX","MIX",post_proc);  CHKERRQ(ierr);
     ierr = post_proc.writeFile(out_file.c_str()); CHKERRQ(ierr);
     PetscFunctionReturn(0);
   }

◆ solveLinearProblem()

PetscErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )

solve problem

Returns: error code

Definition at line 497 of file MixTransportElement.hpp.

                                       {
 
     PetscFunctionBegin;
 
     ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
     ierr = VecZeroEntries(F); CHKERRQ(ierr);
     ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
     ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
     ierr = VecZeroEntries(D0); CHKERRQ(ierr);
     ierr = VecGhostUpdateBegin(D0,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
     ierr = VecGhostUpdateEnd(D0,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
     ierr = VecZeroEntries(D); CHKERRQ(ierr);
     ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
     ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
 
     ierr = mField.query_interface<VecManager>()->setGlobalGhostVector("MIX",COL,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
 
     // Calculate essential boundary conditions
 
     // clear essential bc indices, it could have dofs from other mesh refinement
     bcIndices.clear();
     // clear operator, just in case if some other operators are left on this element
     feTri.getOpPtrVector().clear();
     // set operator to calculate essential boundary conditions
     feTri.getOpPtrVector().push_back(new OpEvaluateBcOnFluxes(*this,"FLUXES"));
     ierr = mField.loop_finite_elements("MIX","MIX_BCFLUX",feTri); CHKERRQ(ierr);
     ierr = VecGhostUpdateBegin(D0,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
     ierr = VecGhostUpdateEnd(D0,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
     ierr = VecAssemblyBegin(D0); CHKERRQ(ierr);
     ierr = VecAssemblyEnd(D0); CHKERRQ(ierr);
 
     // set operators to calculate matrix and right hand side vectors
     feVol.getOpPtrVector().clear();
     feVol.getOpPtrVector().push_back(new OpL2Source(*this,"VALUES",F));
     feVol.getOpPtrVector().push_back(new OpFluxDivergenceAtGaussPts(*this,"FLUXES"));
     feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this,"VALUES"));
     feVol.getOpPtrVector().push_back(new OpDivTauU_HdivL2(*this,"FLUXES","VALUES",F));
     feVol.getOpPtrVector().push_back(new OpTauDotSigma_HdivHdiv(*this,"FLUXES",Aij,F));
     feVol.getOpPtrVector().push_back(new OpVDivSigma_L2Hdiv(*this,"VALUES","FLUXES",Aij,F));
     ierr = mField.loop_finite_elements("MIX","MIX",feVol); CHKERRQ(ierr);
 
     // calculate right hand side for natural boundary conditions
     feTri.getOpPtrVector().clear();
     feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this,"FLUXES",F));
     ierr = mField.loop_finite_elements("MIX","MIX_BCVALUE",feTri); CHKERRQ(ierr);
 
     // assemble matrices
     ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
     ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
     ierr = VecGhostUpdateBegin(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
     ierr = VecGhostUpdateEnd(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
     ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
     ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
 
     {
       double nrm2_F;
       ierr = VecNorm(F,NORM_2,&nrm2_F); CHKERRQ(ierr);
       PetscPrintf(PETSC_COMM_WORLD,"nrm2_F = %6.4e\n",nrm2_F);
     }
 
     // ierr = MatMultAdd(Aij,D0,F,F); CHKERRQ(ierr);
 
     // for ksp solver vector is moved into rhs side
     // for snes it is left ond the left
     ierr = VecScale(F,-1); CHKERRQ(ierr);
 
     IS essential_bc_ids;
     ierr = getDirichletBCIndices(&essential_bc_ids); CHKERRQ(ierr);
     ierr = MatZeroRowsColumnsIS(Aij,essential_bc_ids,1,D0,F); CHKERRQ(ierr);
     ierr = ISDestroy(&essential_bc_ids); CHKERRQ(ierr);
 
     // {
     //   double norm;
     //   ierr = MatNorm(Aij,NORM_FROBENIUS,&norm); CHKERRQ(ierr);
     //   PetscPrintf(PETSC_COMM_WORLD,"mat norm = %6.4e\n",norm);
     // }
 
     {
       double nrm2_F;
       ierr = VecNorm(F,NORM_2,&nrm2_F); CHKERRQ(ierr);
       PetscPrintf(PETSC_COMM_WORLD,"With BC nrm2_F = %6.4e\n",nrm2_F);
     }
 
     //MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
     //MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
     // std::string wait;
     //std::cin >> wait;
 
     // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
     // std::cin >> wait;
 
     // Solve
     KSP solver;
     ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
     ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
     ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
     ierr = KSPSetUp(solver); CHKERRQ(ierr);
     ierr = KSPSolve(solver,F,D); CHKERRQ(ierr);
     ierr = KSPDestroy(&solver); CHKERRQ(ierr);
 
     {
       double nrm2_D;
       ierr = VecNorm(D,NORM_2,&nrm2_D); CHKERRQ(ierr);
       PetscPrintf(PETSC_COMM_WORLD,"nrm2_D = %6.4e\n",nrm2_D);
     }
     ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
     ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
 
 
     // copy data form vector on mesh
     ierr = mField.query_interface<VecManager>()->setGlobalGhostVector("MIX",COL,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
 
     PetscFunctionReturn(0);
   }

Member Data Documentation

◆ Aij

Mat MixTransport::MixTransportElement::Aij

Definition at line 481 of file MixTransportElement.hpp.

◆ bcIndices

set<int> MixTransport::MixTransportElement::bcIndices

Definition at line 93 of file MixTransportElement.hpp.

◆ D

Vec MixTransport::MixTransportElement::D

Definition at line 480 of file MixTransportElement.hpp.

◆ D0

Vec MixTransport::MixTransportElement::D0

Definition at line 480 of file MixTransportElement.hpp.

◆ divergenceAtGaussPts

VectorDouble MixTransport::MixTransportElement::divergenceAtGaussPts

divergence at integration points on element

Definition at line 90 of file MixTransportElement.hpp.

◆ errorMap

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

Definition at line 1480 of file MixTransportElement.hpp.

◆ F

Vec MixTransport::MixTransportElement::F

Definition at line 480 of file MixTransportElement.hpp.

◆ feTri

MyTriFE MixTransport::MixTransportElement::feTri

Instance of surface element.

Definition at line 73 of file MixTransportElement.hpp.

◆ feVol

MyVolumeFE MixTransport::MixTransportElement::feVol

Instance of volume element.

Definition at line 60 of file MixTransportElement.hpp.

◆ fluxesAtGaussPts

MatrixDouble MixTransport::MixTransportElement::fluxesAtGaussPts

fluxes at integration points on element

Definition at line 91 of file MixTransportElement.hpp.

◆ mField

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

Definition at line 46 of file MixTransportElement.hpp.

◆ setOfBlocks

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

maps block set id with appropriate BlockData

Definition at line 199 of file MixTransportElement.hpp.

◆ sumErrorDiv

double MixTransport::MixTransportElement::sumErrorDiv

Definition at line 1482 of file MixTransportElement.hpp.

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1481 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

Definition at line 1483 of file MixTransportElement.hpp.

◆ valuesAtGaussPts

VectorDouble MixTransport::MixTransportElement::valuesAtGaussPts

values at integration points on element

Definition at line 88 of file MixTransportElement.hpp.

◆ valuesGradientAtGaussPts

MatrixDouble MixTransport::MixTransportElement::valuesGradientAtGaussPts

gradients at integration points on element

Definition at line 89 of file MixTransportElement.hpp.

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

users_modules/basic_finite_elements/mix_transport/src/MixTransportElement.hpp

Classes

Public Member Functions

Public Attributes

Detailed Description

Constructor & Destructor Documentation

◆ MixTransportElement()

◆ ~MixTransportElement()

Member Function Documentation

◆ addFields()

◆ addFiniteElements()

◆ buildProblem()

◆ calculateResidual()

◆ createMatrices()

◆ destroyMatrices()

◆ evaluateError()

◆ getBcOnFluxes()

◆ getBcOnValues()

◆ getDirichletBCIndices()

◆ getResistivity()

◆ getSource()

◆ postProc()

◆ solveLinearProblem()

Member Data Documentation

◆ Aij

◆ bcIndices

◆ D

◆ D0

◆ divergenceAtGaussPts

◆ errorMap

◆ F

◆ feTri

◆ feVol

◆ fluxesAtGaussPts

◆ mField

◆ setOfBlocks

◆ sumErrorDiv

◆ sumErrorFlux

◆ sumErrorJump

◆ valuesAtGaussPts

◆ valuesGradientAtGaussPts