MixTransport::UnsaturatedFlowElement Struct Reference

Implementation of operators, problem and finite elements for unsaturated flow. More...

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

Inheritance diagram for MixTransport::UnsaturatedFlowElement:

Collaboration diagram for MixTransport::UnsaturatedFlowElement:

Classes
struct	BcData
	Class storing information about boundary condition. More...
struct	FaceRule
	Set integration rule to boundary elements. More...
struct	MonitorPostProc
struct	OpDivTauU_HdivL2
struct	OpEvaluateInitiallHead
struct	OpIntegrateFluxes
struct	OpPostProcMaterial
struct	OpResidualFlux
	Assemble flux residual. More...
struct	OpResidualMass
struct	OpRhsBcOnValues
	Evaluate boundary condition at the boundary. More...
struct	OpTauDotSigma_HdivHdiv
struct	OpVDivSigma_L2Hdiv
struct	OpVU_L2L2
struct	postProcessVol
	Post proces method for volume element Assemble vectors and matrices and apply essential boundary conditions. More...
struct	preProcessVol
	Pre-peprocessing Set head pressure rate and get inital essential boundary conditions. More...
struct	VolRule
	Set integration rule to volume elements. More...

Public Types
typedef std::map< int, boost::shared_ptr< GenericMaterial > >	MaterialsDoubleMap
typedef map< int, boost::shared_ptr< BcData > >	BcMap

Public Member Functions
	UnsaturatedFlowElement (MoFEM::Interface &m_field)
	~UnsaturatedFlowElement ()
virtual MoFEMErrorCode	getMaterial (const EntityHandle ent, int &block_id) const
	For given element handle get material block Id. More...
MoFEMErrorCode	getBcOnValues (const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
	Get value on boundary. More...
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 More...
MoFEMErrorCode	addFiniteElements (const std::string &fluxes_name, const std::string &values_name)
	add finite elements More...
MoFEMErrorCode	buildProblem (Range zero_flux_ents, BitRefLevel ref_level=BitRefLevel().set(0))
	Build problem. More...
MoFEMErrorCode	setFiniteElements (ForcesAndSourcesCore::RuleHookFun vol_rule=VolRule(), ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule())
	Create finite element instances. More...
MoFEMErrorCode	createMatrices ()
	Create vectors and matrices. More...
MoFEMErrorCode	destroyMatrices ()
	Delete matrices and vector when no longer needed. More...
MoFEMErrorCode	calculateEssentialBc ()
	Calculate boundary conditions for fluxes. More...
MoFEMErrorCode	calculateInitialPc ()
	Calculate inital pressure head distribution. More...
MoFEMErrorCode	solveProblem (bool set_initial_pc=true)
	solve problem More...
Public Member Functions inherited from MixTransport::MixTransportElement
	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...
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
DM	dM
	Discrete manager for unsaturated flow problem. More...
MaterialsDoubleMap	dMatMap
	materials database More...
BcMap	bcValueMap
	Store boundary condition for head capillary pressure. More...
EntityHandle	lastEvalBcValEnt
int	lastEvalBcBlockValId
BcMap	bcFluxMap
EntityHandle	lastEvalBcFluxEnt
int	lastEvalBcBlockFluxId
boost::shared_ptr< ForcesAndSourcesCore >	feFaceBc
	Elemnet to calculate essential bc. More...
boost::shared_ptr< ForcesAndSourcesCore >	feFaceRhs
	Face element apply natural bc. More...
boost::shared_ptr< ForcesAndSourcesCore >	feVolInitialPc
	Calculate inital boundary conditions. More...
boost::shared_ptr< ForcesAndSourcesCore >	feVolRhs
	Assemble residual vector. More...
boost::shared_ptr< ForcesAndSourcesCore >	feVolLhs
	Assemble tangent matrix. More...
boost::shared_ptr< MethodForForceScaling >	scaleMethodFlux
	Method scaling fluxes. More...
boost::shared_ptr< MethodForForceScaling >	scaleMethodValue
	Method scaling values. More...
boost::shared_ptr< FEMethod >	tsMonitor
boost::shared_ptr< VectorDouble >	headRateAtGaussPts
	Vector keeps head rate. More...
std::vector< int >	bcVecIds
VectorDouble	bcVecVals
VectorDouble	vecValsOnBc
Vec	D1
	Vector with inital head capillary pressure. More...
Vec	ghostFlux
	Ghost Vector of integrated fluxes. More...
Public Attributes inherited from MixTransport::MixTransportElement
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

Implementation of operators, problem and finite elements for unsaturated flow.

Examples

unsaturated_transport.cpp.

Definition at line 102 of file UnsaturatedFlow.hpp.

Member Typedef Documentation

BcMap

typedef map<int, boost::shared_ptr<BcData> > MixTransport::UnsaturatedFlowElement::BcMap

Examples

UnsaturatedFlow.hpp.

Definition at line 152 of file UnsaturatedFlow.hpp.

MaterialsDoubleMap

typedef std::map<int, boost::shared_ptr<GenericMaterial> > MixTransport::UnsaturatedFlowElement::MaterialsDoubleMap

Examples

UnsaturatedFlow.hpp.

Definition at line 118 of file UnsaturatedFlow.hpp.

Constructor & Destructor Documentation

UnsaturatedFlowElement()

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

inline

Examples

UnsaturatedFlow.hpp.

Definition at line 106 of file UnsaturatedFlow.hpp.

      : MixTransportElement(m_field), dM(PETSC_NULL), lastEvalBcValEnt(0),
        lastEvalBcBlockValId(-1), lastEvalBcFluxEnt(0),
        lastEvalBcBlockFluxId(-1) {}

~UnsaturatedFlowElement()

MixTransport::UnsaturatedFlowElement::~UnsaturatedFlowElement ( )

inline

Examples

UnsaturatedFlow.hpp.

Definition at line 111 of file UnsaturatedFlow.hpp.

                            {
    if (dM != PETSC_NULL) {
      CHKERR DMDestroy(&dM);
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
    }
  }

Member Function Documentation

addFields()

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

inline

add fields

Examples

UnsaturatedFlow.hpp.

Definition at line 1136 of file UnsaturatedFlow.hpp.

                                            {
    MoFEMFunctionBegin;
    // Fields
    CHKERR mField.add_field(fluxes, HDIV, DEMKOWICZ_JACOBI_BASE, 1);
    CHKERR mField.add_field(values, L2, AINSWORTH_LEGENDRE_BASE, 1);
    CHKERR mField.add_field(values + "_t", L2, AINSWORTH_LEGENDRE_BASE, 1);
    // CHKERR mField.add_field(fluxes+"_residual",L2,AINSWORTH_LEGENDRE_BASE,1);
 
    // meshset consisting all entities in mesh
    EntityHandle root_set = mField.get_moab().get_root_set();
    // add entities to field
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, 4, "SOIL") != 0)
        continue;
      CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
                                              MBTET, fluxes);
      CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
                                              MBTET, values);
      CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
                                              MBTET, values + "_t");
      // CHKERR mField.add_ents_to_field_by_type(
      //   dMatMap[it->getMeshsetId()]->tEts,MBTET,fluxes+"_residual"
      // );
    }
 
    CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
    CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
    CHKERR mField.set_field_order(root_set, MBTET, values, order);
    CHKERR mField.set_field_order(root_set, MBTET, values + "_t", order);
    // CHKERR mField.set_field_order(root_set,MBTET,fluxes+"_residual",order);
    MoFEMFunctionReturn(0);
  }

addFiniteElements()

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

inline

add finite elements

Examples

UnsaturatedFlow.hpp.

Definition at line 1172 of file UnsaturatedFlow.hpp.

                                                                 {
    MoFEMFunctionBegin;
 
    // Define element "MIX". Note that this element will work with fluxes_name
    // and values_name. This reflect bilinear form for the problem
    CHKERR mField.add_finite_element("MIX", MF_ZERO);
    CHKERR mField.modify_finite_element_add_field_row("MIX", fluxes_name);
    CHKERR mField.modify_finite_element_add_field_col("MIX", fluxes_name);
    CHKERR mField.modify_finite_element_add_field_row("MIX", values_name);
    CHKERR mField.modify_finite_element_add_field_col("MIX", values_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX", fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX", values_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX",
                                                       values_name + "_t");
    // CHKERR
    // mField.modify_finite_element_add_field_data("MIX",fluxes_name+"_residual");
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, 4, "SOIL") != 0)
        continue;
      CHKERR mField.add_ents_to_finite_element_by_type(
          dMatMap[it->getMeshsetId()]->tEts, MBTET, "MIX");
    }
 
    // Define element to integrate natural boundary conditions, i.e. set values.
    CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
    CHKERR mField.modify_finite_element_add_field_row("MIX_BCVALUE",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_col("MIX_BCVALUE",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX_BCVALUE",
                                                       fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX_BCVALUE",
                                                       values_name);
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, 4, "HEAD") != 0)
        continue;
      CHKERR mField.add_ents_to_finite_element_by_type(
          bcValueMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCVALUE");
    }
 
    // Define element to apply essential boundary conditions.
    CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
    CHKERR mField.modify_finite_element_add_field_row("MIX_BCFLUX",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_col("MIX_BCFLUX",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX_BCFLUX",
                                                       fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX_BCFLUX",
                                                       values_name);
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, 4, "FLUX") != 0)
        continue;
      CHKERR mField.add_ents_to_finite_element_by_type(
          bcFluxMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCFLUX");
    }
 
    MoFEMFunctionReturn(0);
  }

buildProblem()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::buildProblem ( Range  zero_flux_ents, BitRefLevel  ref_level = BitRefLevel().set(0)  )

inline

Build problem.

Parameters

ref_level

mesh refinement on which mesh problem you like to built.

Returns

error code

Examples

UnsaturatedFlow.hpp.

Definition at line 1241 of file UnsaturatedFlow.hpp.

                                                                            {
    MoFEMFunctionBegin;
 
    // Build fields
    CHKERR mField.build_fields();
    // Build finite elements
    CHKERR mField.build_finite_elements("MIX");
    CHKERR mField.build_finite_elements("MIX_BCFLUX");
    CHKERR mField.build_finite_elements("MIX_BCVALUE");
    // Build adjacencies of degrees of freedom and elements
    CHKERR mField.build_adjacencies(ref_level);
 
    //  create DM instance
    CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
    // setting that DM is type of DMMOFEM, i.e. MOFEM implementation manages DM
    CHKERR DMSetType(dM, "DMMOFEM");
    // mesh is portioned, each process keeps only part of problem
    CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
    // creates problem in DM
    CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
    // discretised problem creates square matrix (that makes some optimizations)
    CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
    // set DM options from command line
    CHKERR DMSetFromOptions(dM);
    // add finite elements
    CHKERR DMMoFEMAddElement(dM, "MIX");
    CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
    CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
    // constructor data structures
    CHKERR DMSetUp(dM);
 
    // remove zero flux dofs
    CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
        "MIX", "FLUXES", zero_flux_ents);
 
    PetscSection section;
    CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
    CHKERR DMSetSection(dM, section);
    // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
    CHKERR PetscSectionDestroy(&section);
 
    MoFEMFunctionReturn(0);
  }

calculateEssentialBc()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::calculateEssentialBc ( )

inline

Calculate boundary conditions for fluxes.

Returns

Error code

Definition at line 1650 of file UnsaturatedFlow.hpp.

                                        {
    MoFEMFunctionBegin;
    // clear vectors
    CHKERR VecZeroEntries(D0);
    CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
    // clear essential bc indices, it could have dofs from other mesh refinement
    bcIndices.clear();
    // set operator to calculate essential boundary conditions
    CHKERR DMoFEMLoopFiniteElements(dM, "MIX_BCFLUX", feFaceBc);
    CHKERR VecAssemblyBegin(D0);
    CHKERR VecAssemblyEnd(D0);
    CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
    double norm2D0;
    CHKERR VecNorm(D0, NORM_2, &norm2D0);
    // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
    PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
    MoFEMFunctionReturn(0);
  }

calculateInitialPc()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::calculateInitialPc ( )

inline

Calculate inital pressure head distribution.

Returns

Error code

Definition at line 1675 of file UnsaturatedFlow.hpp.

                                      {
    MoFEMFunctionBegin;
    // clear vectors
    CHKERR VecZeroEntries(D1);
    CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
    // Calculate initial pressure head on each element
    CHKERR DMoFEMLoopFiniteElements(dM, "MIX", feVolInitialPc);
    // Assemble vector
    CHKERR VecAssemblyBegin(D1);
    CHKERR VecAssemblyEnd(D1);
    CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
    // Calculate norm
    double norm2D1;
    CHKERR VecNorm(D1, NORM_2, &norm2D1);
    // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
    PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
    MoFEMFunctionReturn(0);
  }

createMatrices()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::createMatrices ( )

inline

Create vectors and matrices.

Returns

Error code

Definition at line 1616 of file UnsaturatedFlow.hpp.

                                  {
    MoFEMFunctionBegin;
    CHKERR DMCreateMatrix(dM, &Aij);
    CHKERR DMCreateGlobalVector(dM, &D0);
    CHKERR VecDuplicate(D0, &D1);
    CHKERR VecDuplicate(D0, &D);
    CHKERR VecDuplicate(D0, &F);
    int ghosts[] = {0};
    int nb_locals = mField.get_comm_rank() == 0 ? 1 : 0;
    int nb_ghosts = mField.get_comm_rank() > 0 ? 1 : 0;
    CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
                          &ghostFlux);
    MoFEMFunctionReturn(0);
  }

destroyMatrices()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::destroyMatrices ( )

inline

Delete matrices and vector when no longer needed.

Returns

error code

Definition at line 1635 of file UnsaturatedFlow.hpp.

                                   {
    MoFEMFunctionBegin;
    CHKERR MatDestroy(&Aij);
    CHKERR VecDestroy(&D);
    CHKERR VecDestroy(&D0);
    CHKERR VecDestroy(&D1);
    CHKERR VecDestroy(&F);
    CHKERR VecDestroy(&ghostFlux);
    MoFEMFunctionReturn(0);
  }

getBcOnFluxes()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::getBcOnFluxes ( const EntityHandle  ent, const double  x, const double  y, const double  z, double &  flux  )

inlinevirtual

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 from MixTransport::MixTransportElement.

Examples

UnsaturatedFlow.hpp.

Definition at line 210 of file UnsaturatedFlow.hpp.

                                                                             {
    MoFEMFunctionBegin;
    int block_id = -1;
    if (lastEvalBcFluxEnt == ent) {
      block_id = lastEvalBcBlockFluxId;
    } else {
      for (BcMap::iterator it = bcFluxMap.begin(); it != bcFluxMap.end();
           it++) {
        if (it->second->eNts.find(ent) != it->second->eNts.end()) {
          block_id = it->first;
        }
      }
      lastEvalBcFluxEnt = ent;
      lastEvalBcBlockFluxId = block_id;
    }
    if (block_id >= 0) {
      if (bcFluxMap.at(block_id)->hookFun) {
        flux = bcFluxMap.at(block_id)->hookFun(x, y, z);
      } else {
        flux = bcFluxMap.at(block_id)->fixValue;
      }
    } else {
      flux = 0;
    }
    MoFEMFunctionReturn(0);
  }

getBcOnValues()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::getBcOnValues ( const EntityHandle  ent, const int  gg, const double  x, const double  y, const double  z, double &  value  )

inlinevirtual

Get value on boundary.

Parameters

ent	entity handle
gg	number of integration point
x	x-coordinate
y	y-coordinate
z	z-coordinate
value	returned value

Returns

error code

Reimplemented from MixTransport::MixTransportElement.

Examples

UnsaturatedFlow.hpp.

Definition at line 168 of file UnsaturatedFlow.hpp.

                                              {
    MoFEMFunctionBegin;
    int block_id = -1;
    if (lastEvalBcValEnt == ent) {
      block_id = lastEvalBcBlockValId;
    } else {
      for (BcMap::iterator it = bcValueMap.begin(); it != bcValueMap.end();
           it++) {
        if (it->second->eNts.find(ent) != it->second->eNts.end()) {
          block_id = it->first;
        }
      }
      lastEvalBcValEnt = ent;
      lastEvalBcBlockValId = block_id;
    }
    if (block_id >= 0) {
      if (bcValueMap.at(block_id)->hookFun) {
        value = bcValueMap.at(block_id)->hookFun(x, y, z);
      } else {
        value = bcValueMap.at(block_id)->fixValue;
      }
    } else {
      value = 0;
    }
    MoFEMFunctionReturn(0);
  }

getMaterial()

virtual MoFEMErrorCode MixTransport::UnsaturatedFlowElement::getMaterial ( const EntityHandle  ent, int &  block_id  ) const

inlinevirtual

For given element handle get material block Id.

Parameters

ent	finite element entity handle
block_id	reference to returned block id

Returns

error code

Examples

UnsaturatedFlow.hpp.

Definition at line 127 of file UnsaturatedFlow.hpp.

                                                          {
    MoFEMFunctionBegin;
    for (MaterialsDoubleMap::const_iterator mit = dMatMap.begin();
         mit != dMatMap.end(); mit++) {
      if (mit->second->tEts.find(ent) != mit->second->tEts.end()) {
        block_id = mit->first;
        MoFEMFunctionReturnHot(0);
      }
    }
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Element not found, no material data");
    MoFEMFunctionReturn(0);
  }

setFiniteElements()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::setFiniteElements ( ForcesAndSourcesCore::RuleHookFun  vol_rule = VolRule(), ForcesAndSourcesCore::RuleHookFun  face_rule = FaceRule()  )

inline

Create finite element instances.

Parameters

vol_rule	integration rule for volume element
face_rule	integration rule for boundary element

Returns

error code

Examples

UnsaturatedFlow.hpp.

Definition at line 1446 of file UnsaturatedFlow.hpp.

                                                                            {
    MoFEMFunctionBegin;
 
    // create finite element instances
    feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
        new FaceElementForcesAndSourcesCore(mField));
    feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
        new FaceElementForcesAndSourcesCore(mField));
    feVolInitialPc = boost::shared_ptr<ForcesAndSourcesCore>(
        new VolumeElementForcesAndSourcesCore(mField));
    feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
        new VolumeElementForcesAndSourcesCore(mField));
    feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
        new VolumeElementForcesAndSourcesCore(mField));
    // set integration rule to elements
    feFaceBc->getRuleHook = face_rule;
    feFaceRhs->getRuleHook = face_rule;
    feVolInitialPc->getRuleHook = vol_rule;
    feVolLhs->getRuleHook = vol_rule;
    feVolRhs->getRuleHook = vol_rule;
    // set function hook for finite element postprocessing stage
    feVolRhs->preProcessHook = preProcessVol(*this, feVolRhs);
    feVolLhs->preProcessHook = preProcessVol(*this, feVolLhs);
    feVolRhs->postProcessHook = postProcessVol(*this, feVolRhs);
    feVolLhs->postProcessHook = postProcessVol(*this, feVolLhs);
 
    // Set Piola transform
    CHKERR AddHOOps<2, 3, 3>::add(feFaceBc->getOpPtrVector(), {HDIV});
    CHKERR AddHOOps<2, 3, 3>::add(feFaceRhs->getOpPtrVector(), {HDIV});
 
    // create method for setting history for fluxes on boundary
    scaleMethodFlux = boost::shared_ptr<MethodForForceScaling>(
        new TimeForceScale("-flux_history", false));
 
    // create method for setting history for presure heads on boundary
    scaleMethodValue = boost::shared_ptr<MethodForForceScaling>(
        new TimeForceScale("-head_history", false));
 
    // Set operator to calculate essential boundary conditions
    feFaceBc->getOpPtrVector().push_back(
        new OpEvaluateBcOnFluxes(*this, "FLUXES"));
 
    // Set operator to calculate initial capillary pressure
    CHKERR AddHOOps<3, 3, 3>::add(feVolInitialPc->getOpPtrVector(), {HDIV, L2});
    feVolInitialPc->getOpPtrVector().push_back(
        new OpEvaluateInitiallHead(*this, "VALUES"));
 
    // set residual face from Neumann terms, i.e. applied pressure
    feFaceRhs->getOpPtrVector().push_back(
        new OpRhsBcOnValues(*this, "FLUXES", scaleMethodValue));
    // set residual finite element operators
    headRateAtGaussPts = boost::make_shared<VectorDouble>();
 
 
    // resAtGaussPts = boost::make_shared<VectorDouble>();
    CHKERR AddHOOps<3, 3, 3>::add(feVolRhs->getOpPtrVector(), {HDIV, L2});
    feVolRhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
        string("VALUES") + "_t", headRateAtGaussPts, MBTET));
    feVolRhs->getOpPtrVector().push_back(
        new OpValuesAtGaussPts(*this, "VALUES"));
    feVolRhs->getOpPtrVector().push_back(
        new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
    feVolRhs->getOpPtrVector().push_back(new OpResidualFlux(*this, "FLUXES"));
    feVolRhs->getOpPtrVector().push_back(new OpResidualMass(*this, "VALUES"));
    feVolRhs->getOpPtrVector().back().opType =
        ForcesAndSourcesCore::UserDataOperator::OPROW;
    // set tangent matrix finite element operators
    CHKERR AddHOOps<3, 3, 3>::add(feVolLhs->getOpPtrVector(), {HDIV, L2});
    feVolLhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
        string("VALUES") + "_t", headRateAtGaussPts, MBTET));
    // feVolLhs->getOpPtrVector().push_back(
    //   new
    //   OpCalculateScalarFieldValues(string("FLUXES")+"_residual",resAtGaussPts,MBTET)
    // );
    feVolLhs->getOpPtrVector().push_back(
        new OpValuesAtGaussPts(*this, "VALUES"));
    feVolLhs->getOpPtrVector().push_back(
        new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
    feVolLhs->getOpPtrVector().push_back(
        new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
    feVolLhs->getOpPtrVector().push_back(new OpVU_L2L2(*this, "VALUES"));
    feVolLhs->getOpPtrVector().push_back(
        new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
    feVolLhs->getOpPtrVector().push_back(
        new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
 
    // Adding finite elements to DM, time solver will ask for it to assemble
    // tangent matrices and residuals
    boost::shared_ptr<FEMethod> null;
    CHKERR DMMoFEMTSSetIFunction(dM, "MIX_BCVALUE", feFaceRhs, null, null);
    CHKERR DMMoFEMTSSetIFunction(dM, "MIX", feVolRhs, null, null);
    CHKERR DMMoFEMTSSetIJacobian(dM, "MIX", feVolLhs, null, null);
 
    // setting up post-processing
 
 
    auto get_post_process_ele = [&]() {
      auto post_process = boost::make_shared<PostProcEle>(mField);
 
      CHKERR AddHOOps<3, 3, 3>::add(post_process->getOpPtrVector(), {HDIV, L2});
 
      auto values_ptr = boost::make_shared<VectorDouble>();
      auto grad_ptr = boost::make_shared<MatrixDouble>();
      auto flux_ptr = boost::make_shared<MatrixDouble>();
 
      post_process->getOpPtrVector().push_back(
          new OpCalculateScalarFieldValues("VALUES", values_ptr));
      post_process->getOpPtrVector().push_back(
          new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
      post_process->getOpPtrVector().push_back(
          new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
 
      using OpPPMap = OpPostProcMapInMoab<3, 3>;
 
      post_process->getOpPtrVector().push_back(
 
          new OpPPMap(post_process->getPostProcMesh(),
                      post_process->getMapGaussPts(),
 
                      {{"VALUES", values_ptr}},
 
                      {{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
 
                      {}, {}
 
                      ));
 
      return post_process;
    };
 
    auto post_process = get_post_process_ele();
 
    post_process->getOpPtrVector().push_back(
        new OpValuesAtGaussPts(*this, "VALUES"));
    post_process->getOpPtrVector().push_back(
        new OpPostProcMaterial(*this, post_process->getPostProcMesh(),
                               post_process->getMapGaussPts(), "VALUES"));
 
    // Integrate fluxes on boundary
    boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
    flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
        new FaceElementForcesAndSourcesCore(mField));
    CHKERR AddHOOps<2, 3, 3>::add(flux_integrate->getOpPtrVector(), {HDIV});
    flux_integrate->getOpPtrVector().push_back(
        new OpIntegrateFluxes(*this, "FLUXES"));
    int frequency = 1;
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Monitor post proc", "none");
    CHKERR PetscOptionsInt(
        "-how_often_output",
        "frequency how often results are dumped on hard disk", "", frequency,
        &frequency, NULL);
    ierr = PetscOptionsEnd();
    CHKERRG(ierr);
 
    tsMonitor = boost::shared_ptr<FEMethod>(
        new MonitorPostProc(*this, post_process, flux_integrate, frequency));
    TsCtx *ts_ctx;
    DMMoFEMGetTsCtx(dM, &ts_ctx);
    ts_ctx->getPostProcessMonitor().push_back(tsMonitor);
    MoFEMFunctionReturn(0);
  }

solveProblem()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::solveProblem ( bool  set_initial_pc = true )

inline

solve problem

Returns

error code

Definition at line 1700 of file UnsaturatedFlow.hpp.

                                                          {
    MoFEMFunctionBegin;
    if (set_initial_pc) {
      // Set initial head
      CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
    }
 
    // Initiate vector from data on the mesh
    CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_FORWARD);
 
    // Create time solver
    TS ts;
    CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
    // Use backward Euler method
    CHKERR TSSetType(ts, TSBEULER);
    // Set final time
    double ftime = 1;
    CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
    // Setup solver from commabd line
    CHKERR TSSetFromOptions(ts);
    // Set DM to TS
    CHKERR TSSetDM(ts, dM);
#if PETSC_VERSION_GE(3, 7, 0)
    CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
#endif
    // Set-up monitor
    TsCtx *ts_ctx;
    DMMoFEMGetTsCtx(dM, &ts_ctx);
    CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
 
    // This add SNES monitor, to show error by fields. It is dirty trick
    // to add monitor, so code is hidden from doxygen
    CHKERR TSSetSolution(ts, D);
    CHKERR TSSetUp(ts);
    SNES snes;
    CHKERR TSGetSNES(ts, &snes);
 
#if PETSC_VERSION_GE(3, 7, 0)
    {
      PetscViewerAndFormat *vf;
      CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
                                        PETSC_VIEWER_DEFAULT, &vf);
      CHKERR SNESMonitorSet(
          snes,
          (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
                             void *))SNESMonitorFields,
          vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
    }
#else
    {
      CHKERR SNESMonitorSet(snes,
                            (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
                                               void *))SNESMonitorFields,
                            0, 0);
    }
#endif
 
    CHKERR TSSolve(ts, D);
 
    // Get statistic form TS and print it
    CHKERR TSGetTime(ts, &ftime);
    PetscInt steps, snesfails, rejects, nonlinits, linits;
    CHKERR TSGetTimeStepNumber(ts, &steps);
    CHKERR TSGetSNESFailures(ts, &snesfails);
    CHKERR TSGetStepRejections(ts, &rejects);
    CHKERR TSGetSNESIterations(ts, &nonlinits);
    CHKERR TSGetKSPIterations(ts, &linits);
    PetscPrintf(PETSC_COMM_WORLD,
                "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
                "%D, linits %D\n",
                steps, rejects, snesfails, ftime, nonlinits, linits);
 
    MoFEMFunctionReturn(0);
  }

Member Data Documentation

bcFluxMap

BcMap MixTransport::UnsaturatedFlowElement::bcFluxMap

Examples

UnsaturatedFlow.hpp.

Definition at line 197 of file UnsaturatedFlow.hpp.

bcValueMap

BcMap MixTransport::UnsaturatedFlowElement::bcValueMap

Store boundary condition for head capillary pressure.

Examples

UnsaturatedFlow.hpp.

Definition at line 153 of file UnsaturatedFlow.hpp.

bcVecIds

std::vector<int> MixTransport::UnsaturatedFlowElement::bcVecIds

Examples

UnsaturatedFlow.hpp.

Definition at line 1323 of file UnsaturatedFlow.hpp.

bcVecVals

VectorDouble MixTransport::UnsaturatedFlowElement::bcVecVals

Examples

UnsaturatedFlow.hpp.

Definition at line 1324 of file UnsaturatedFlow.hpp.

D1

Vec MixTransport::UnsaturatedFlowElement::D1

Vector with inital head capillary pressure.

Examples

UnsaturatedFlow.hpp.

Definition at line 1609 of file UnsaturatedFlow.hpp.

dM

DM MixTransport::UnsaturatedFlowElement::dM

Discrete manager for unsaturated flow problem.

Examples

UnsaturatedFlow.hpp.

Definition at line 104 of file UnsaturatedFlow.hpp.

dMatMap

MaterialsDoubleMap MixTransport::UnsaturatedFlowElement::dMatMap

materials database

Examples

UnsaturatedFlow.hpp.

Definition at line 119 of file UnsaturatedFlow.hpp.

feFaceBc

boost::shared_ptr<ForcesAndSourcesCore> MixTransport::UnsaturatedFlowElement::feFaceBc

Elemnet to calculate essential bc.

Examples

UnsaturatedFlow.hpp.

Definition at line 1287 of file UnsaturatedFlow.hpp.

feFaceRhs

boost::shared_ptr<ForcesAndSourcesCore> MixTransport::UnsaturatedFlowElement::feFaceRhs

Face element apply natural bc.

Examples

UnsaturatedFlow.hpp.

Definition at line 1289 of file UnsaturatedFlow.hpp.

feVolInitialPc

boost::shared_ptr<ForcesAndSourcesCore> MixTransport::UnsaturatedFlowElement::feVolInitialPc

Calculate inital boundary conditions.

Examples

UnsaturatedFlow.hpp.

Definition at line 1291 of file UnsaturatedFlow.hpp.

feVolLhs

boost::shared_ptr<ForcesAndSourcesCore> MixTransport::UnsaturatedFlowElement::feVolLhs

Assemble tangent matrix.

Examples

UnsaturatedFlow.hpp.

Definition at line 1294 of file UnsaturatedFlow.hpp.

feVolRhs

boost::shared_ptr<ForcesAndSourcesCore> MixTransport::UnsaturatedFlowElement::feVolRhs

Assemble residual vector.

Examples

UnsaturatedFlow.hpp.

Definition at line 1293 of file UnsaturatedFlow.hpp.

ghostFlux

Vec MixTransport::UnsaturatedFlowElement::ghostFlux

Ghost Vector of integrated fluxes.

Examples

UnsaturatedFlow.hpp.

Definition at line 1610 of file UnsaturatedFlow.hpp.

headRateAtGaussPts

boost::shared_ptr<VectorDouble> MixTransport::UnsaturatedFlowElement::headRateAtGaussPts

Vector keeps head rate.

Examples

UnsaturatedFlow.hpp.

Definition at line 1303 of file UnsaturatedFlow.hpp.

lastEvalBcBlockFluxId

int MixTransport::UnsaturatedFlowElement::lastEvalBcBlockFluxId

Examples

UnsaturatedFlow.hpp.

Definition at line 199 of file UnsaturatedFlow.hpp.

lastEvalBcBlockValId

int MixTransport::UnsaturatedFlowElement::lastEvalBcBlockValId

Examples

UnsaturatedFlow.hpp.

Definition at line 156 of file UnsaturatedFlow.hpp.

lastEvalBcFluxEnt

EntityHandle MixTransport::UnsaturatedFlowElement::lastEvalBcFluxEnt

Examples

UnsaturatedFlow.hpp.

Definition at line 198 of file UnsaturatedFlow.hpp.

lastEvalBcValEnt

EntityHandle MixTransport::UnsaturatedFlowElement::lastEvalBcValEnt

Examples

UnsaturatedFlow.hpp.

Definition at line 155 of file UnsaturatedFlow.hpp.

scaleMethodFlux

boost::shared_ptr<MethodForForceScaling> MixTransport::UnsaturatedFlowElement::scaleMethodFlux

Method scaling fluxes.

Examples

UnsaturatedFlow.hpp.

Definition at line 1296 of file UnsaturatedFlow.hpp.

scaleMethodValue

boost::shared_ptr<MethodForForceScaling> MixTransport::UnsaturatedFlowElement::scaleMethodValue

Method scaling values.

Examples

UnsaturatedFlow.hpp.

Definition at line 1298 of file UnsaturatedFlow.hpp.

tsMonitor

boost::shared_ptr<FEMethod> MixTransport::UnsaturatedFlowElement::tsMonitor

Element used by TS monitor to postprocess results at time step

Examples

UnsaturatedFlow.hpp.

Definition at line 1299 of file UnsaturatedFlow.hpp.

vecValsOnBc

VectorDouble MixTransport::UnsaturatedFlowElement::vecValsOnBc

Examples

UnsaturatedFlow.hpp.

Definition at line 1324 of file UnsaturatedFlow.hpp.

---

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

• tutorials/cor-0to1/src/UnsaturatedFlow.hpp