Mix transport problem. More...

#include "tutorials/cor-0to1/src/MixTransportElement.hpp"

Inheritance diagram for MixTransport::MixTransportElement:

Collaboration diagram for MixTransport::MixTransportElement:

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

virtual	~MixTransportElement ()
	destructor

MoFEMErrorCode	getDirichletBCIndices (IS *is)
	get dof indices where essential boundary conditions are applied

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

virtual MoFEMErrorCode	getResistivity (const EntityHandle ent, const double x, const double y, const double z, MatrixDouble3by3 &inv_k)
	natural (Dirichlet) boundary conditions (set values)

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

virtual MoFEMErrorCode	getBcOnFluxes (const EntityHandle ent, const double x, const double y, const double z, double &flux)
	essential (Neumann) boundary condition (set fluxes)

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

MoFEMErrorCode	addFiniteElements (const std::string &fluxes_name, const std::string &values_name)
	add finite elements

MoFEMErrorCode	buildProblem (BitRefLevel &ref_level)
	Build problem.

MoFEMErrorCode	postProc (const string out_file)
	Post process results.

MoFEMErrorCode	createMatrices ()
	create matrices

MoFEMErrorCode	solveLinearProblem ()
	solve problem

MoFEMErrorCode	calculateResidual ()
	calculate residual

MoFEMErrorCode	evaluateError ()
	Calculate error on elements.

MoFEMErrorCode	destroyMatrices ()
	destroy matrices

Public Attributes
MoFEM::Interface &	mField

MyVolumeFE	feVol
	Instance of volume element.

MyTriFE	feTri
	Instance of surface element.

VectorDouble	valuesAtGaussPts
	values at integration points on element

MatrixDouble	valuesGradientAtGaussPts
	gradients at integration points on element

VectorDouble	divergenceAtGaussPts
	divergence at integration points on element

MatrixDouble	fluxesAtGaussPts
	fluxes at integration points on element

set< int >	bcIndices

std::map< int, BlockData >	setOfBlocks
	maps block set id with appropriate BlockData

Vec	D

Vec	D0

Vec	F

Mat	Aij

map< double, EntityHandle >	errorMap

double	sumErrorFlux

double	sumErrorDiv

double	sumErrorJump

Detailed Description

Mix transport problem.

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

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

General problem have form,

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

Definition at line 34 of file MixTransportElement.hpp.

Constructor & Destructor Documentation

◆ MixTransportElement()

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

inline

construction of this data structure

Examples: UnsaturatedFlow.hpp.

Definition at line 70 of file MixTransportElement.hpp.

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

MixTransport::MixTransportElement::mField

MoFEM::Interface & mField

Definition MixTransportElement.hpp:36

MixTransport::MixTransportElement::feTri

MyTriFE feTri

Instance of surface element.

Definition MixTransportElement.hpp:65

MixTransport::MixTransportElement::feVol

MyVolumeFE feVol

Instance of volume element.

Definition MixTransportElement.hpp:51

◆ ~MixTransportElement()

virtual MixTransport::MixTransportElement::~MixTransportElement ( )

inlinevirtual

destructor

Definition at line 76 of file MixTransportElement.hpp.

76{}

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

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

Returns: error code

Definition at line 194 of file MixTransportElement.hpp.

                                            {
    MoFEMFunctionBegin;
    // Fields
    CHKERR mField.add_field(fluxes, HDIV, DEMKOWICZ_JACOBI_BASE, 1);
    CHKERR mField.add_field(values, L2, AINSWORTH_LEGENDRE_BASE, 1);
 
    // meshset consisting all entities in mesh
    EntityHandle root_set = mField.get_moab().get_root_set();
    // add entities to field
    CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, fluxes);
    CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, values);
    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);
    MoFEMFunctionReturn(0);
  }

◆ addFiniteElements()

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

inline

add finite elements

Definition at line 213 of file MixTransportElement.hpp.

                                                                 {
    MoFEMFunctionBegin;
 
    // 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
    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);
 
    // Define finite element to integrate over skeleton, we need that to
    // evaluate error
    CHKERR mField.add_finite_element("MIX_SKELETON", MF_ZERO);
    CHKERR mField.modify_finite_element_add_field_row("MIX_SKELETON",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_col("MIX_SKELETON",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX_SKELETON",
                                                       fluxes_name);
 
    // 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)) {
 
      Mat_Thermal temp_data;
      CHKERR it->getAttributeDataStructure(temp_data);
      setOfBlocks[it->getMeshsetId()].cOnductivity =
          temp_data.data.Conductivity;
      setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
      CHKERR mField.get_moab().get_entities_by_type(
          it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
      CHKERR mField.add_ents_to_finite_element_by_type(
          setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
 
      Range skeleton;
      CHKERR mField.get_moab().get_adjacencies(
          setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
          moab::Interface::UNION);
      CHKERR mField.add_ents_to_finite_element_by_type(skeleton, MBTRI,
                                                       "MIX_SKELETON");
    }
 
    // 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);
 
    // 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);
    MoFEMFunctionReturn(0);
  }

◆ buildProblem()

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

inline

Build problem.

Parameters

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

Returns: error code

Definition at line 296 of file MixTransportElement.hpp.

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

◆ calculateResidual()

MoFEMErrorCode MixTransport::MixTransportElement::calculateResidual ( )

inline

calculate residual

Definition at line 617 of file MixTransportElement.hpp.

                                     {
    MoFEMFunctionBegin;
    CHKERR VecZeroEntries(F);
    CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecAssemblyBegin(F);
    CHKERR VecAssemblyEnd(F);
    // calculate residuals
    feVol.getOpPtrVector().clear();
    CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
    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_NULLPTR, F));
    feVol.getOpPtrVector().push_back(
        new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", PETSC_NULLPTR, F));
    CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
    feTri.getOpPtrVector().clear();
    feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
    CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
    CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(F);
    CHKERR VecAssemblyEnd(F);
    // CHKERR VecAXPY(F,-1.,D0);
    // CHKERR MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULLPTR,F);
    {
      std::vector<int> ids;
      ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
      std::vector<double> vals(ids.size(), 0);
      CHKERR VecSetValues(F, ids.size(), &*ids.begin(), &*vals.begin(),
                          INSERT_VALUES);
      CHKERR VecAssemblyBegin(F);
      CHKERR VecAssemblyEnd(F);
    }
    {
      double nrm2_F;
      CHKERR VecNorm(F, NORM_2, &nrm2_F);
      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");
      }
    }
    MoFEMFunctionReturn(0);
  }

◆ createMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::createMatrices ( )

inline

create matrices

Definition at line 479 of file MixTransportElement.hpp.

                                  {
    MoFEMFunctionBegin;
    CHKERR mField.getInterface<MatrixManager>()
        ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("MIX", &Aij);
    CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
    CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
    CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
    MoFEMFunctionReturn(0);
  }

◆ destroyMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::destroyMatrices ( )

inline

destroy matrices

Definition at line 705 of file MixTransportElement.hpp.

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

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

                                 {
    MoFEMFunctionBegin;
    errorMap.clear();
    sumErrorFlux = 0;
    sumErrorDiv = 0;
    sumErrorJump = 0;
    feTri.getOpPtrVector().clear();
    feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
    CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
                                       mField.get_comm_size());
    feVol.getOpPtrVector().clear();
    CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
    feVol.getOpPtrVector().push_back(
        new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
    feVol.getOpPtrVector().push_back(
        new OpValuesGradientAtGaussPts(*this, "VALUES"));
    feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
    CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
                                       mField.get_comm_size());
    const Problem *problem_ptr;
    CHKERR mField.get_problem("MIX", &problem_ptr);
    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);
    MoFEMFunctionReturn(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

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

                                                     {
    MoFEMFunctionBeginHot;
    flux = 0;
    MoFEMFunctionReturnHot(0);
  }

◆ 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

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

                                                                      {
    MoFEMFunctionBeginHot;
    value = 0;
    MoFEMFunctionReturnHot(0);
  }

◆ getDirichletBCIndices()

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

inline

get dof indices where essential boundary conditions are applied

Parameters

is indices

Returns: error code

Definition at line 92 of file MixTransportElement.hpp.

                                               {
    MoFEMFunctionBegin;
    std::vector<int> ids;
    ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
    IS is_local;
    CHKERR ISCreateGeneral(mField.get_comm(), ids.size(),
                           ids.empty() ? PETSC_NULLPTR : &ids[0],
                           PETSC_COPY_VALUES, &is_local);
    CHKERR ISAllGather(is_local, is);
    CHKERR ISDestroy(&is_local);
    MoFEMFunctionReturn(0);
  }

◆ 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

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

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

◆ 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

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

                                                 {
    MoFEMFunctionBeginHot;
    flux = 0;
    MoFEMFunctionReturnHot(0);
  }

◆ postProc()

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

inline

Post process results.

Returns: error code

Definition at line 433 of file MixTransportElement.hpp.

                                                 {
    MoFEMFunctionBegin;
 
    PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
        mField);
 
    CHKERR AddHOOps<3, 3, 3>::add(post_proc.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_proc.getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("VALUES", values_ptr));
    post_proc.getOpPtrVector().push_back(
        new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
    post_proc.getOpPtrVector().push_back(
        new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
 
    using OpPPMap = OpPostProcMapInMoab<3, 3>;
 
    post_proc.getOpPtrVector().push_back(
 
        new OpPPMap(post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
 
                    {{"VALUES", values_ptr}},
 
                    {{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
 
                    {}, {}
 
                    ));
 
    post_proc.getOpPtrVector().push_back(new OpPostProc(
        post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
 
    CHKERR mField.loop_finite_elements("MIX", "MIX", post_proc);
 
    CHKERR post_proc.writeFile(out_file.c_str());
    MoFEMFunctionReturn(0);
  }

◆ solveLinearProblem()

MoFEMErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )

inline

solve problem

Returns: error code

Definition at line 493 of file MixTransportElement.hpp.

                                      {
    MoFEMFunctionBegin;
    CHKERR MatZeroEntries(Aij);
    CHKERR VecZeroEntries(F);
    CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecZeroEntries(D0);
    CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecZeroEntries(D);
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
 
    CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
        "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
 
    // 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();
    CHKERR AddHOOps<2, 3, 3>::add(feTri.getOpPtrVector(), {HDIV});
    // set operator to calculate essential boundary conditions
    feTri.getOpPtrVector().push_back(new OpEvaluateBcOnFluxes(*this, "FLUXES"));
    CHKERR mField.loop_finite_elements("MIX", "MIX_BCFLUX", feTri);
    CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(D0);
    CHKERR VecAssemblyEnd(D0);
 
    // set operators to calculate matrix and right hand side vectors
    feVol.getOpPtrVector().clear();
    CHKERR AddHOOps<3, 3, 3>::add(feVol.getOpPtrVector(), {HDIV, L2});
    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));
    CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
 
    // calculate right hand side for natural boundary conditions
    feTri.getOpPtrVector().clear();
    CHKERR AddHOOps<2, 3, 3>::add(feTri.getOpPtrVector(), {HDIV});
    feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
    CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
 
    // assemble matrices
    CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
    CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(F);
    CHKERR VecAssemblyEnd(F);
 
    {
      double nrm2_F;
      CHKERR VecNorm(F, NORM_2, &nrm2_F);
      PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
    }
 
    // CHKERR MatMultAdd(Aij,D0,F,F); ;
 
    // for ksp solver vector is moved into rhs side
    // for snes it is left ond the left
    CHKERR VecScale(F, -1);
 
    IS essential_bc_ids;
    CHKERR getDirichletBCIndices(&essential_bc_ids);
    CHKERR MatZeroRowsColumnsIS(Aij, essential_bc_ids, 1, D0, F);
    CHKERR ISDestroy(&essential_bc_ids);
 
    // {
    //   double norm;
    //   CHKERR MatNorm(Aij,NORM_FROBENIUS,&norm); ;
    //   PetscPrintf(PETSC_COMM_WORLD,"mat norm = %6.4e\n",norm);
    // }
 
    {
      double nrm2_F;
      CHKERR VecNorm(F, NORM_2, &nrm2_F);
      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;
    CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
    CHKERR KSPSetOperators(solver, Aij, Aij);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetUp(solver);
    CHKERR KSPSolve(solver, F, D);
    CHKERR KSPDestroy(&solver);
 
    {
      double nrm2_D;
      CHKERR VecNorm(D, NORM_2, &nrm2_D);
      ;
      PetscPrintf(PETSC_COMM_WORLD, "nrm2_D = %6.4e\n", nrm2_D);
    }
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
 
    // copy data form vector on mesh
    CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
        "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
 
    MoFEMFunctionReturn(0);
  }

Member Data Documentation

◆ Aij

Mat MixTransport::MixTransportElement::Aij

Definition at line 476 of file MixTransportElement.hpp.

◆ bcIndices

set<int> MixTransport::MixTransportElement::bcIndices

Examples: UnsaturatedFlow.hpp.

Definition at line 85 of file MixTransportElement.hpp.

◆ D

Vec MixTransport::MixTransportElement::D

Definition at line 475 of file MixTransportElement.hpp.

◆ D0

Vec MixTransport::MixTransportElement::D0

Examples: UnsaturatedFlow.hpp.

Definition at line 475 of file MixTransportElement.hpp.

◆ divergenceAtGaussPts

VectorDouble MixTransport::MixTransportElement::divergenceAtGaussPts

divergence at integration points on element

Examples: UnsaturatedFlow.hpp.

Definition at line 82 of file MixTransportElement.hpp.

◆ errorMap

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

Definition at line 1365 of file MixTransportElement.hpp.

◆ F

Vec MixTransport::MixTransportElement::F

Definition at line 475 of file MixTransportElement.hpp.

◆ feTri

MyTriFE MixTransport::MixTransportElement::feTri

Instance of surface element.

Definition at line 65 of file MixTransportElement.hpp.

◆ feVol

MyVolumeFE MixTransport::MixTransportElement::feVol

Instance of volume element.

Definition at line 51 of file MixTransportElement.hpp.

◆ fluxesAtGaussPts

MatrixDouble MixTransport::MixTransportElement::fluxesAtGaussPts

fluxes at integration points on element

Examples: UnsaturatedFlow.hpp.

Definition at line 83 of file MixTransportElement.hpp.

◆ mField

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

Examples: UnsaturatedFlow.hpp.

Definition at line 36 of file MixTransportElement.hpp.

◆ setOfBlocks

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

maps block set id with appropriate BlockData

Definition at line 185 of file MixTransportElement.hpp.

◆ sumErrorDiv

double MixTransport::MixTransportElement::sumErrorDiv

Definition at line 1367 of file MixTransportElement.hpp.

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1366 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

Definition at line 1368 of file MixTransportElement.hpp.

◆ valuesAtGaussPts

VectorDouble MixTransport::MixTransportElement::valuesAtGaussPts

values at integration points on element

Examples: UnsaturatedFlow.hpp.

Definition at line 78 of file MixTransportElement.hpp.

◆ valuesGradientAtGaussPts

MatrixDouble MixTransport::MixTransportElement::valuesGradientAtGaussPts

gradients at integration points on element

Definition at line 80 of file MixTransportElement.hpp.

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

tutorials/cor-0to1/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