v0.14.0
COR-10: Navier-Stokes equation

# Introduction

Navier-Stokes equations (NSE), governing the motion of a viscous fluid, are used in various applications: from simulations of the flow in blood vessels to studies of the air flow around aeroplane wings and rotor blades, scaling up to models of ocean and atmospheric currents. Even in the case of an incompressible steady flow, NSE are non-linear due to the effect of the inertia, which is more pronounced in case of a higher Reynolds number. In this tutorial we discuss the implementation of a viscous fluid model in MoFEM using hierarchical basis functions. This approach permits us to locally increase the order of approximation, enforcing conformity across finite element boundaries, without the need to change the implementation of an element. Moreover, the requirement of different approximation orders for primal (velocity) and dual (pressure) variables, necessary for a simulation of the flow using the mixed formulation, can be easily satisfied.

# Problem statement

An incompressible isoviscous steady-state flow in a domain $$\Omega$$ is governed by the following equations:

\begin{align} \label{eq:balance_momentum} \rho \left(\mathbf{u}\cdot\nabla\right)\mathbf{u} - \mu\nabla^2\mathbf{u} + \nabla p &= \mathbf{f},\\ \label{eq:cont} \nabla\cdot\mathbf{u} &= 0, \end{align}

where \eqref{eq:balance_momentum} is the set of Navier-Stokes equations, representing the balance of the momentum, and \eqref{eq:cont} is the continuity equation; $$\mathbf{u}=[u_x, u_y, u_z]^\intercal$$ is the velocity field, $$p$$ is the hydrostatic pressure field, $$\rho$$ is the fluid mass density, $$\mu$$ is fluid viscosity and $\mathbf{f}$ is the density of external forces. The boundary value problem complements the above equations by the Dirichlet and Neumann conditions on the boundary $$\partial\Omega$$:

\begin{align} \label{eq:bc_d}\mathbf{u} = \mathbf{u}_D & \;\text{on}\; \Gamma_D,\\ \label{eq:bc_n}\mathbf{n}\cdot \left[-p\mathbf{I}+\mu\left(\nabla\mathbf{u} + \nabla\mathbf{u}^\intercal\right)\right] = \mathbf{g}_N & \;\text{on} \;\Gamma_N, \end{align}

where $$\mathbf{u}_D$$ is the prescribed velocity on the part of the boundary $$\Gamma_D\subset\partial\Omega$$, and $$\mathbf{g}_N$$ is the prescribed traction vector on the part of the boundary $$\Gamma_N\subset\partial\Omega$$, $$\mathbf{n}$$ is an outward normal.

In this tutorial we will solve the problem of the fluid flow around a rigid sphere (see Figure 1) of a radius $$r$$ positioned in the centre of a cubic domain, the side length $$2l$$ of which is considered sufficiently large compared to $$r$$, so that a uniform far-field velocity on the exterior boundaries is valid (note that the body forces are neglected). Exploiting the symmetry of the problem, we will use a quarter of the domain in our simulations.

Figure 1: (a) Finite-element mesh used for simulation of fluid flow around a rigid sphere. (b) Sketch of the problem set-up on a section $z=0$ of the mesh.

Note that due to the structure of differential operators in the Navier-Stokes equations \eqref{eq:balance_momentum}, the fluid pressure has to be specified on a part of boundary in order to obtain a unique solution. To achieve that, we will combine Dirichlet and Neumann condition on the outlet boundary of the domain, see see Figure 1. In the case when the flow is perpendicular to the boundary ( $$u_y = u_z = 0$$), which is enforced by the Dirichlet boundary condition, the hydrostatic fluid pressure equals to the normal traction on such boundary, see Theorem 1 in [11].

# Non-dimensionalization and scaling

Reynolds number is introduced for the problems governed by the Navier-Stokes equations as a measure of the ratio of inertial forces to viscous forces:

$$$\mathcal{R} = \frac{\rho U L }{\mu},$$$

where $$U$$ is the scale for the velocity and $$L$$ is a relevant length scale. For the problem of the fluid flow around a sphere, considered in this tutorial, the far-field velocity plays the role of the velocity scale, while $$L=2r$$, where $$r$$ is the radius of the sphere.

When viscous forces dominate over the inertia forces, $$\mathcal{R} \ll 1$$, the non-linear term in \eqref{eq:balance_momentum} can be neglected, simplifying NSE down to Stokes equations. However, often this is not the case, i.e. the nonlinear term may become dominant over the viscous forces. In these cases non-dimensionalization (scaling) of Navier-Stokes equations is very useful, which permits to decrease the number of coefficients to only one – Reynolds number $$\mathcal{R}$$. The following scales can be used for the case of flow with relatively low Reynolds number (e.g. less than 100).

Non-dimensionalization of Navier-Stokes equations
Physical quantity Scale Dimensionless variable
Length $$L$$ $$\mathbf{x}^{*}=\mathbf{x}\:/\:L, \quad \nabla^{*}(\cdot) = L\nabla(\cdot)$$
Velocity $$U$$ $$\mathbf{u}^{*}=\mathbf{u}\:/\:U$$
Pressure $$P=\frac{\mu U}{L}$$ $$p^{*}=p\:/\:P$$

Mote that the scale for pressure depends on the scales for the length and the velocity. Using these scales, Navier-stokes equations (1) - (2) can be written in the following dimensionless form (considering zero external forces):

\begin{align} \label{eq:balance_momentum_dim_less} \mathcal{R} \left(\mathbf{u}^{*}\cdot\nabla^{*}\right)\mathbf{u}^{*} - (\nabla^{*})^2\mathbf{u}^{*} + \nabla^{*} p^{*} &= 0,\\ \label{eq:cont_dim_less} \nabla^{*}\cdot\mathbf{u}^{*} &= 0, \end{align}

Note that now the influence of the non-linearity (inertia terms) on the problem is fully controlled by the value of the Reynolds number. On the one hand, if $$\mathcal{R} \ll 1$$, or if one simply wants to consider the Stokes equation, the nonlinear term can be dropped. On the other hand, if the Reynolds number is sufficiently high, then the non-linearity can become strong and may pose problems for convergence (see below discussion of the linearisation of the problem for the Newton-Raphson method). In this case a certain iterative procedure can be used to obtain a solution. Indeed, to find the solution for a given Reynolds number, a number of intermediate problems can be solved for range of smaller $$\mathcal{R}$$, eventually reaching the original value. Such technique will also be used in this tutorial. Note also that below we will drop the * symbol, implying that during computation all considered quantities are dimensionless, and a transformation to the dimensional variables is performed before post-processing of the results takes place.

# Finite-element formulation

## Weak form

The weak statement of the problem \eqref{eq:balance_momentum}-\eqref{eq:bc_n} reads: Find a vector field $$\mathbf{u} \in \mathbf{H}^1(\Omega)$$ and a scalar field $$p \in L^2(\Omega)$$, such that for any test functions $$\mathbf{v} \in \mathbf{H}^1_0(\Omega)$$ and $$q \in L^2(\Omega)$$:

\begin{align} \label{eq:weak} \mathcal{R}\int\limits_\Omega \left(\mathbf{u}\cdot\nabla\right)\mathbf{u} \cdot\mathbf{v} \,d\Omega + \int\limits_{\Omega}\nabla\mathbf{u}\mathbin{:}\nabla\mathbf{v}\, d\Omega - \int\limits_{\Omega}p\, \nabla\cdot\mathbf{v} \, d\Omega &= \int\limits_{\Omega}\mathbf{f}\cdot\mathbf{v}\,d\Omega + \int\limits_{\Gamma_N}\mathbf{g}_N\cdot\mathbf{v}\,d\Gamma_N,\\ - \int\limits_{\Omega}q\, \nabla\cdot\mathbf{u} \, d\Omega &= 0 \end{align}

Upon finite-element discretization of the domain $$\Omega$$, we consider interpolation of both unknown fields introducing shape functions on each element: $$\label{eq:shape} u_i = \sum\limits_{\alpha=1}^{n_{\mathbf{u}}} N_{\alpha}\, u_i^\alpha, \quad p = \sum\limits_{\beta=1}^{n_p} \Phi_{\beta}\, p^\beta; \quad v_i = \sum\limits_{\alpha=1}^{n_{\mathbf{u}}} N_{\alpha}\, v_i^\alpha, \quad q = \sum\limits_{\beta=1}^{n_p} \Phi_{\beta}\, q^\beta,$$ where $$n_{\mathbf{u}}$$ is the number of shape functions associated with the velocity field, and $n_p$ is the similar number for the pressure field. Using the hierarchical basis approximation, the vector of the shape functions can be decomposed into four sub-vectors, consisting of shape functions associated with element's entities: vertices, edges, faces and the volume of the element, e.g. for the velocity field:

$$$\label{eq:shape_hier} \mathbf{N}^\textit{el} = \left[N_1, \ldots N_\alpha, \ldots N_{n_{\mathbf{u}}}\right]^\intercal = \left[\mathbf{N}^\textit{ver}, \mathbf{N}^\textit{edge}, \mathbf{N}^\textit{face}, \mathbf{N}^\textit{vol}\right]^\intercal.$$$

## Note on the choice of spaces for the mixed problem

Note that in the weak statement velocity and pressure fields belong to different spaces: $$H^1$$ and $$L^2$$, respectively. Two different approaches can be used for implementing such mixed formulation in the finite-element framework. One can use generalised Taylor-Hood elements, which feature continuous ( $$H^1$$ ) approximation of both velocity and pressure fields, however, in order to enforce stability, the approximations functions for the pressure field should be one order lower than those for the velocity field. Furthermore, Taylor-Hood elements impose certain constraints on the mesh, see [77]. Alternatively, a discontinuous approximation for pressure can be used, however, in this case the difference between the orders of approximation functions for velocity and for pressure should be 2.

# Implementation

The example class has necessary fields to store input parameters and internal data structures, as well as a set of functions for setting up and solving the problem:

commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
}
CHKERR SNESDestroy(&snes);
CHKERR DMDestroy(&dM);
}
private:
PetscBool isPartitioned = PETSC_FALSE;
int orderVelocity = 2; // default approximation orderVelocity
int orderPressure = 1; // default approximation orderPressure
int numHoLevels = 0;
PetscBool isDiscontPressure = PETSC_FALSE;
PetscBool isHoGeometry = PETSC_TRUE;
double pressureScale = 1.0;
double lengthScale = 1.0;
double velocityScale = 1.0;
double reNumber = 1.0;
double density;
double viscosity;
PetscBool isStokesFlow = PETSC_FALSE;
int numSteps = 1;
double lambdaStep = 1.0;
double lambda = 0.0;
int step;
DM dM;
SNES snes;
boost::shared_ptr<NavierStokesElement::CommonData> commonData;
SmartPetscObj<Vec> D;
SmartPetscObj<Vec> F;
SmartPetscObj<Mat> Aij;
boost::shared_ptr<VolumeElementForcesAndSourcesCore> feLhsPtr;
boost::shared_ptr<VolumeElementForcesAndSourcesCore> feRhsPtr;
boost::shared_ptr<FaceElementForcesAndSourcesCore> feDragPtr;
boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feDragSidePtr;
boost::shared_ptr<DirichletDisplacementBc> dirichletBcPtr;
boost::ptr_map<std::string, NeumannForcesSurface> neumannForces;
boost::shared_ptr<PostProcVol> fePostProcPtr;
boost::shared_ptr<PostProcFace> fePostProcDragPtr;
};

The set of functions used in this tutorials in the following:

The function ExampleNavierStokes::runProblem is executed from the main function:

int main(int argc, char *argv[]) {
const char param_file[] = "param_file.petsc";
// Initialise MoFEM
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
try {
// Create mesh database
moab::Core mb_instance; // create database
moab::Interface &moab = mb_instance; // create interface to database
// Create MoFEM database and link it to MoAB
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
NavierStokesExample ex(m_field);
CHKERR ex.runProblem();
}
// finish work cleaning memory, getting statistics, etc
return 0;
}

## Reading mesh and input parameters

The workflow of solving a problem using the finite-element method starts from reading the mesh and other input parameters:

char mesh_file_name[255];
PetscBool flg_mesh_file;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "NAVIER_STOKES problem",
"none");
CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
mesh_file_name, 255, &flg_mesh_file);
CHKERR PetscOptionsBool("-my_is_partitioned", "is partitioned?", "",
isPartitioned, &isPartitioned, PETSC_NULL);
CHKERR PetscOptionsInt("-my_order_u", "approximation orderVelocity", "",
orderVelocity, &orderVelocity, PETSC_NULL);
CHKERR PetscOptionsInt("-my_order_p", "approximation orderPressure", "",
orderPressure, &orderPressure, PETSC_NULL);
CHKERR PetscOptionsInt("-my_num_ho_levels",
"number of higher order boundary levels", "",
numHoLevels, &numHoLevels, PETSC_NULL);
CHKERR PetscOptionsBool("-my_discont_pressure", "discontinuous pressure", "",
CHKERR PetscOptionsBool("-my_ho_geometry", "use second order geometry", "",
isHoGeometry, &isHoGeometry, PETSC_NULL);
CHKERR PetscOptionsScalar("-my_length_scale", "length scale", "", lengthScale,
&lengthScale, PETSC_NULL);
CHKERR PetscOptionsScalar("-my_velocity_scale", "velocity scale", "",
velocityScale, &velocityScale, PETSC_NULL);
CHKERR PetscOptionsBool("-my_stokes_flow", "stokes flow", "", isStokesFlow,
&isStokesFlow, PETSC_NULL);
CHKERR PetscOptionsInt("-my_step_num", "number of steps", "", numSteps,
&numSteps, PETSC_NULL);
ierr = PetscOptionsEnd();
if (flg_mesh_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
// Read whole mesh or part of it if partitioned
if (isPartitioned == PETSC_TRUE) {
// This is a case of distributed mesh and algebra. In that case each
// processor keeps only part of the problem.
const char *option;
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
CHKERR mField.get_moab().load_file(mesh_file_name, 0, option);
} else {
// In this case, we have distributed algebra, i.e. assembly of vectors and
// matrices is in parallel, but whole mesh is stored on all processors.
// snes and matrix scales well, however problem set-up of problem is
// not fully parallel.
const char *option;
option = "";
CHKERR mField.get_moab().load_file(mesh_file_name, 0, option);
}
bitLevel.set(0);
CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(0, 3,
}

Once the mesh is read, we find the set of tetrahedra corresponding to the computational domain. Furthermore, we identify the fluid-structure interface, i.e. the set of triangles corresponding to the surface of the sphere:

if (bit->getName().compare(0, 9, "MAT_FLUID") == 0) {
const int id = bit->getMeshsetId();
CHKERR mField.get_moab().get_entities_by_dimension(
bit->getMeshset(), 3, commonData->setOfBlocksData[id].eNts, true);
std::vector<double> attributes;
bit->getAttributes(attributes);
if (attributes.size() < 2) {
SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"should be at least 2 attributes but is %d",
attributes.size());
}
commonData->setOfBlocksData[id].iD = id;
commonData->setOfBlocksData[id].fluidViscosity = attributes[0];
commonData->setOfBlocksData[id].fluidDensity = attributes[1];
viscosity = commonData->setOfBlocksData[id].fluidViscosity;
density = commonData->setOfBlocksData[id].fluidDensity;
}
}
if (bit->getName().compare(0, 9, "INT_SOLID") == 0) {
Range tets, tet;
const int id = bit->getMeshsetId();
CHKERR mField.get_moab().get_entities_by_type(
bit->getMeshset(), MBTRI, commonData->setOfFacesData[id].eNts, true);
solidFaces.merge(commonData->setOfFacesData[id].eNts);
commonData->setOfFacesData[id].eNts, 3, true, tets,
moab::Interface::UNION);
tet = Range(tets.front(), tets.front());
for (auto &bit : commonData->setOfBlocksData) {
if (bit.second.eNts.contains(tet)) {
commonData->setOfFacesData[id].fluidViscosity =
bit.second.fluidViscosity;
commonData->setOfFacesData[id].fluidDensity = bit.second.fluidDensity;
commonData->setOfFacesData[id].iD = id;
break;
}
}
if (commonData->setOfFacesData[id].fluidViscosity < 0) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Cannot find a fluid block adjacent to a given solid face");
}
}
}
}

## Setting up fields, finite elements and the problem itself

Once the domain and the fluid-structure interface are identified, we can setup the problem. First, we define the velocity and pressure fields. Note that for velocity we use space H1, while for pressure the user can choose between using also H1 space (continuous approximation, Taylor-Hood element), or a discontinuous approximation using L2 space. Note that in the former case, the approximation functions for pressure has to be one order lower than that of the velocity (e.g. order 2 for velocity and order 1 for pressure), while the the latter case the difference between the two has to be 2 orders (e.g. order 3 for velocity and order 1 for pressure).

Furthermore, the properties of hierarchical basis functions can be used to locally increase the approximation order for the elements on the fluid-structure interface, where the gradients of the velocity can be expected to be the highest due to the no-slip boundary condition on the surface of the sphere.

} else {
}
3);
CHKERR mField.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
CHKERR mField.set_field_order(0, MBVERTEX, "VELOCITY", 1);
CHKERR mField.set_field_order(0, MBEDGE, "VELOCITY", orderVelocity);
CHKERR mField.set_field_order(0, MBTRI, "VELOCITY", orderVelocity);
CHKERR mField.set_field_order(0, MBTET, "VELOCITY", orderVelocity);
CHKERR mField.set_field_order(0, MBVERTEX, "PRESSURE", 1);
CHKERR mField.set_field_order(0, MBEDGE, "PRESSURE", orderPressure);
CHKERR mField.set_field_order(0, MBTRI, "PRESSURE", orderPressure);
}
CHKERR mField.set_field_order(0, MBTET, "PRESSURE", orderPressure);
if (numHoLevels > 0) {
std::vector<Range> levels(numHoLevels);
Range ents;
ents.merge(solidFaces);
for (int ll = 0; ll != numHoLevels; ll++) {
Range verts;
CHKERR mField.get_moab().get_connectivity(ents, verts, true);
for (auto d : {1, 2, 3}) {
CHKERR mField.get_moab().get_adjacencies(verts, d, false, ents,
moab::Interface::UNION);
}
levels[ll] = subtract(ents, ents.subset_by_type(MBVERTEX));
}
for (int ll = numHoLevels - 1; ll >= 1; ll--) {
levels[ll] = subtract(levels[ll], levels[ll - 1]);
}
int add_order = 1;
for (int ll = numHoLevels - 1; ll >= 0; ll--) {
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
levels[ll]);
CHKERR mField.set_field_order(levels[ll], "VELOCITY",
CHKERR mField.set_field_order(levels[ll], "PRESSURE",
else
CHKERR mField.set_field_order(levels[ll].subset_by_type(MBTET),
"PRESSURE", orderPressure + add_order);
}
}
CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
// Set 2nd order of approximation of geometry
if (isHoGeometry) {
Range ents, edges;
CHKERR mField.get_moab().get_entities_by_dimension(0, 3, ents);
CHKERR mField.get_moab().get_adjacencies(ents, 1, false, edges,
moab::Interface::UNION);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(edges);
CHKERR mField.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
}
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
"VELOCITY");
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
"PRESSURE");
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
"MESH_NODE_POSITIONS");
}
Projection10NodeCoordsOnField ent_method_material(mField,
"MESH_NODE_POSITIONS");
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
}

After that the finite elements operating on the set above fields can be defined. In particular, we define elements for solving Navier-Stokes equations, computation of the drag traction and force on the surface of the surface, and, finally, for computing the contribution of Neumann (natural) boundary conditions.

// Add finite element (this defines element, declaration comes later)
CHKERR NavierStokesElement::addElement(mField, "NAVIER_STOKES", "VELOCITY",
"PRESSURE", "MESH_NODE_POSITIONS");
CHKERR NavierStokesElement::addElement(mField, "DRAG", "VELOCITY", "PRESSURE",
"MESH_NODE_POSITIONS", 2, &solidFaces);
// build finite elements
// build adjacencies between elements and degrees of freedom
}

Once the finite elements are defined, we will need to create element instances and push necessary operators to their pipelines. First, we setup an element for solving Navier-Stokes (or if one prefers, Stokes) equations. Note that right-hand side and left-hand side elements are considered separately. In case of Navier-Stokes equations, operators are pushed to elements pipelines using function NavierStokesElement::setNavierStokesOperators. In particular, the following operators are pushed to right-hand side:

while for the left-hand side we push these operators:

Upon that, we setup in the same way elements for computing the drag traction and the drag force acting on the sphere using operators NavierStokesElement::OpCalcDragTraction and NavierStokesElement::OpCalcDragForce, respectively. Additionally, operators for the Neumann boundary conditions are set using the function MetaNeumannForces::setMomentumFluxOperators and the element for the Dirichlet boundary conditions is created as well. Finally, operators for post-processing output are created: in the volume and on the fluid-solid interface.

feLhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
feRhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feLhsPtr, true, false, false,
true);
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feRhsPtr, true, false, false,
true);
feDragPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
feDragSidePtr = boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>(
new VolumeElementForcesAndSourcesCoreOnSide(mField));
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feDragSidePtr, true, false, false,
true);
if (isStokesFlow) {
feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
} else {
feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
}
"NAVIER_STOKES", "VELOCITY",
"PRESSURE", commonData);
dirichletBcPtr = boost::shared_ptr<DirichletDisplacementBc>(
new DirichletDisplacementBc(mField, "VELOCITY", Aij, D, F));
// dirichletBcPtr->snes_ctx = FEMethod::CTX_SNESNONE;
dirichletBcPtr->snes_x = D;
// Assemble pressure and traction forces
NULL, "VELOCITY");
// for postprocessing:
fePostProcPtr = boost::make_shared<PostProcVol>(mField);
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fePostProcPtr, true, false, false,
true);
auto v_ptr = boost::make_shared<MatrixDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
auto pos_ptr = boost::make_shared<MatrixDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
fePostProcPtr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("VELOCITY", v_ptr));
fePostProcPtr->getOpPtrVector().push_back(
fePostProcPtr->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PRESSURE", p_ptr));
fePostProcPtr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
using OpPPMap = OpPostProcMapInMoab<3, 3>;
fePostProcPtr->getOpPtrVector().push_back(
new OpPPMap(
fePostProcPtr->getPostProcMesh(), fePostProcPtr->getMapGaussPts(),
{{"PRESSURE", p_ptr}},
{{"VELOCITY", v_ptr}, {"MESH_NODE_POSITIONS", pos_ptr}},
{}
)
);
fePostProcDragPtr = boost::make_shared<PostProcFace>(mField);
fePostProcDragPtr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
fePostProcDragPtr->getOpPtrVector().push_back(
new OpPPMap(
fePostProcDragPtr->getPostProcMesh(),
fePostProcDragPtr->getMapGaussPts(),
{},
{{"MESH_NODE_POSITIONS", pos_ptr}},
{},
{}
)
);
fePostProcDragPtr, feDragSidePtr, "NAVIER_STOKES", "VELOCITY", "PRESSURE",
commonData);
}

Note that the class for enforcing Dirichlet boundary conditions DirichletDisplacementBc used in the snippet above requires the algebraic data structures for storing components of the right-hand side vector, left-hand side matrix and the solution vector. Therefore, these data structures need to be created beforehand using PETSc:

CHKERR VecZeroEntries(F);
CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecZeroEntries(D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
CHKERR MatZeroEntries(Aij);
}

However, one may note that creation of PETSCc vectors and matrices requires the instance of the MoFEM Discrete Manager (DM), which needs a certain setup as well:

DMType dm_name = "DM_NAVIER_STOKES";
// Register DM problem
CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
CHKERR DMSetType(dM, dm_name);
// Create DM instance
// Configure DM form line command options s
CHKERR DMSetFromOptions(dM);
// Add elements to dM (only one here)
if (mField.check_finite_element("PRESSURE_FE"))
// setup the DM
CHKERR DMSetUp(dM);
}

Finally, the last preparatory step required is the setup of the SNES functionality:

boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
neumannForces.begin();
for (; mit != neumannForces.end(); mit++) {
CHKERR DMMoFEMSNESSetFunction(dM, mit->first.c_str(),
&mit->second->getLoopFe(), NULL, NULL);
}
boost::shared_ptr<FEMethod> null_fe;
CHKERR DMMoFEMSNESSetFunction(dM, "NAVIER_STOKES", feRhsPtr, null_fe,
null_fe);
null_fe);
CHKERR DMMoFEMSNESSetJacobian(dM, "NAVIER_STOKES", feLhsPtr, null_fe,
null_fe);
SnesCtx *snes_ctx;
// create snes nonlinear solver
CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
CHKERR SNESSetFromOptions(snes);
}

## Solving the problem and post-processing the results

auto scale_problem = [&](double U, double L, double P) {
CHKERR mField.getInterface<FieldBlas>()->fieldScale(L,
"MESH_NODE_POSITIONS");
ProjectionFieldOn10NodeTet ent_method_on_10nodeTet(
mField, "MESH_NODE_POSITIONS", true, true);
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_on_10nodeTet);
ent_method_on_10nodeTet.setNodes = false;
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_on_10nodeTet);
CHKERR mField.getInterface<FieldBlas>()->fieldScale(U, "VELOCITY");
CHKERR mField.getInterface<FieldBlas>()->fieldScale(P, "PRESSURE");
};
CHKERR scale_problem(1.0 / velocityScale, 1.0 / lengthScale,
1.0 / pressureScale);
step = 0;
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Viscosity: %6.4e\n", viscosity);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Density: %6.4e\n", density);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Length scale: %6.4e\n", lengthScale);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Velocity scale: %6.4e\n",
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Pressure scale: %6.4e\n",
if (isStokesFlow) {
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Re number: 0 (Stokes flow)\n");
} else {
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Re number: %6.4e\n",
}
while (lambda < 1.0 - 1e-12) {
if (isStokesFlow) {
reNumber = 0.0;
for (auto &bit : commonData->setOfBlocksData) {
bit.second.inertiaCoef = 0.0;
bit.second.viscousCoef = 1.0;
}
} else {
for (auto &bit : commonData->setOfBlocksData) {
bit.second.inertiaCoef = reNumber;
bit.second.viscousCoef = 1.0;
}
}
CHKERR PetscPrintf(
PETSC_COMM_WORLD,
"Step: %d | Lambda: %6.4e | Inc: %6.4e | Re number: %6.4e \n", step,
CHKERR VecAssemblyBegin(D);
CHKERR VecAssemblyEnd(D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR SNESSolve(snes, PETSC_NULL, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToGlobalVector(dM, D, INSERT_VALUES, SCATTER_REVERSE);
CHKERR scale_problem(1.0 / velocityScale, 1.0 / lengthScale,
1.0 / pressureScale);
step++;
}
}
string out_file_name;
{
std::ostringstream stm;
stm << "out_" << step << ".h5m";
out_file_name = stm.str();
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n",
out_file_name.c_str());
CHKERR fePostProcPtr->writeFile(out_file_name.c_str());
}
CHKERR VecZeroEntries(commonData->pressureDragForceVec);
CHKERR VecZeroEntries(commonData->shearDragForceVec);
CHKERR VecZeroEntries(commonData->totalDragForceVec);
auto get_vec_data = [&](auto vec, std::array<double, 3> &data) {
CHKERR VecAssemblyBegin(vec);
CHKERR VecAssemblyEnd(vec);
const double *array;
if (mField.get_comm_rank() == 0) {
for (int i : {0, 1, 2})
data[i] = array[i];
}
};
std::array<double, 3> pressure_drag;
std::array<double, 3> shear_drag;
std::array<double, 3> total_drag;
CHKERR get_vec_data(commonData->pressureDragForceVec, pressure_drag);
CHKERR get_vec_data(commonData->shearDragForceVec, shear_drag);
CHKERR get_vec_data(commonData->totalDragForceVec, total_drag);
if (mField.get_comm_rank() == 0) {
MOFEM_LOG_C("SELF", Sev::inform,
"Pressure drag force: [ %6.4e, %6.4e, %6.4e ]",
pressure_drag[0], pressure_drag[1], pressure_drag[2]);
MOFEM_LOG_C("SELF", Sev::inform,
"Shear drag force: [ %6.4e, %6.4e, %6.4e ]", shear_drag[0],
shear_drag[1], shear_drag[2]);
MOFEM_LOG_C("SELF", Sev::inform,
"Total drag force: [ %6.4e, %6.4e, %6.4e ]", total_drag[0],
total_drag[1], total_drag[2]);
}
{
std::ostringstream stm;
stm << "out_drag_" << step << ".h5m";
out_file_name = stm.str();
CHKERR PetscPrintf(PETSC_COMM_WORLD, "out file %s\n",
out_file_name.c_str());
}
}

# Running the example

We will consider two cases: first, flow governed by Stokes equation, i.e. neglecting the nonlinear terms, and then we will solve full set Navier-Stokes equations. For both cases we will use the same mesh. We will start by partitioning this mesh:

mofem_part -my_file sphere_test.cub -my_nparts 4 -output_file sphere_test_4.h5m

while mofem_part can be found in \$HOME/mofem_install/um/build_release/tools/. Here we partitioned the mesh in 4 parts in order to be executed on 4 cores in a distributed-memory parallel environment. Note that any number of parts less or equal to number of physically available cores can be used.

## Stokes equation

mpirun -np 4 ./navier_stokes -my_file sphere_test_4.h5m \
-my_order_u 2 \
-my_order_p 1 \
-my_discont_pressure 0 \
-my_num_ho_levels 1 \
-my_ho_geometry 1 \
-my_step_num 1 \
-my_velocity_scale 0.1 \
-my_length_scale 2 \
-my_stokes_flow 1 \
-my_is_partitioned 1

## Navier-Stokes equations

mpirun -np 4 ./navier_stokes -my_file sphere_test_4.h5m \
-my_order_u 2 \
-my_order_p 1 \
-my_discont_pressure 0 \
-my_num_ho_levels 1 \
-my_ho_geometry 1 \
-my_step_num 2 \
-my_velocity_scale 0.1 \
-my_length_scale 2 \
-my_stokes_flow 0 \
-my_is_partitioned 1
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
NavierStokesExample::density
double density
Definition: navier_stokes.cpp:49
NavierStokesExample::NavierStokesExample
NavierStokesExample(MoFEM::Interface &m_field)
Definition: navier_stokes.cpp:24
NavierStokesExample::runProblem
MoFEMErrorCode runProblem()
[Example Navier Stokes]
Definition: navier_stokes.cpp:96
NavierStokesExample::F
SmartPetscObj< Vec > F
Definition: navier_stokes.cpp:67
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
Definition: SurfacePressure.cpp:1974
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
NavierStokesExample::feDragSidePtr
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > feDragSidePtr
Definition: navier_stokes.cpp:74
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
sdf_hertz.d
float d
Definition: sdf_hertz.py:5
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:764
NavierStokesExample::feDragPtr
boost::shared_ptr< FaceElementForcesAndSourcesCore > feDragPtr
Definition: navier_stokes.cpp:73
NavierStokesExample::solidFaces
Range solidFaces
Definition: navier_stokes.cpp:58
NavierStokesExample
[Example Navier Stokes]
Definition: navier_stokes.cpp:22
NavierStokesExample::dirichletBcPtr
boost::shared_ptr< DirichletDisplacementBc > dirichletBcPtr
Definition: navier_stokes.cpp:76
NavierStokesElement::setPostProcDragOperators
static MoFEMErrorCode setPostProcDragOperators(boost::shared_ptr< PostProcFace > postProcDragPtr, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for post processing output of drag traction.
Definition: NavierStokesElement.cpp:184
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
NavierStokesExample::isDiscontPressure
PetscBool isDiscontPressure
Definition: navier_stokes.cpp:42
NavierStokesExample::Aij
SmartPetscObj< Mat > Aij
Definition: navier_stokes.cpp:68
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
NavierStokesExample::mField
MoFEM::Interface & mField
Definition: navier_stokes.cpp:35
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
out_file_name
char out_file_name[255]
Definition: initial_diffusion.cpp:53
NavierStokesExample::fePostProcDragPtr
boost::shared_ptr< PostProcFace > fePostProcDragPtr
Definition: navier_stokes.cpp:80
NavierStokesExample::findBlocks
MoFEMErrorCode findBlocks()
Definition: navier_stokes.cpp:182
NavierStokesExample::dM
DM dM
Definition: navier_stokes.cpp:61
DirichletDisplacementBc
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
NavierStokesExample::orderVelocity
int orderVelocity
Definition: navier_stokes.cpp:39
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1003
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:501
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
NavierStokesElement::VolRule
Set integration rule to volume elements.
Definition: NavierStokesElement.hpp:263
NavierStokesExample::isPartitioned
PetscBool isPartitioned
Definition: navier_stokes.cpp:37
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::CoreInterface::rebuild_database
virtual MoFEMErrorCode rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
NavierStokesExample::~NavierStokesExample
~NavierStokesExample()
Definition: navier_stokes.cpp:27
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
NavierStokesExample::defineFiniteElements
MoFEMErrorCode defineFiniteElements()
[Setup fields]
Definition: navier_stokes.cpp:334
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1018
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
NavierStokesExample::step
int step
Definition: navier_stokes.cpp:56
NavierStokesExample::orderPressure
int orderPressure
Definition: navier_stokes.cpp:40
MetaNeumannForces::setMomentumFluxOperators
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Definition: SurfacePressure.cpp:2069
NavierStokesExample::viscosity
double viscosity
Definition: navier_stokes.cpp:50
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
DM_NO_ELEMENT
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
NavierStokesExample::setupElementInstances
MoFEMErrorCode setupElementInstances()
[Setup algebraic structures]
Definition: navier_stokes.cpp:404
NavierStokesExample::lambda
double lambda
Definition: navier_stokes.cpp:55
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:193
NavierStokesExample::feLhsPtr
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhsPtr
Definition: navier_stokes.cpp:70
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:539
NavierStokesExample::setupFields
MoFEMErrorCode setupFields()
[Find blocks]
Definition: navier_stokes.cpp:239
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1098
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
NavierStokesExample::bitLevel
BitRefLevel bitLevel
Definition: navier_stokes.cpp:59
NavierStokesExample::D
SmartPetscObj< Vec > D
Definition: navier_stokes.cpp:66
MoFEM::DMMoFEMSNESSetJacobian
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMoFEM.cpp:763
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:550
NavierStokesExample::lengthScale
double lengthScale
Definition: navier_stokes.cpp:46
static double lambda
Definition: NavierStokesElement.hpp:104
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
NavierStokesExample::neumannForces
boost::ptr_map< std::string, NeumannForcesSurface > neumannForces
Definition: navier_stokes.cpp:77
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
[Run problem]
Definition: navier_stokes.cpp:112
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
FaceElementForcesAndSourcesCore
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
NavierStokesExample::pressureScale
double pressureScale
Definition: navier_stokes.cpp:45
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const string element_name, const string velocity_field_name, const string pressure_field_name, const string mesh_field_name, const int dim=3, Range *ents=nullptr)
Setting up elements.
Definition: NavierStokesElement.hpp:129
NavierStokesElement::setCalcDragOperators
static MoFEMErrorCode setCalcDragOperators(boost::shared_ptr< FaceElementForcesAndSourcesCore > dragFe, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for calculating drag force on the solid surface.
Definition: NavierStokesElement.cpp:150
NavierStokesExample::reNumber
double reNumber
Definition: navier_stokes.cpp:48
NavierStokesExample::snes
SNES snes
Definition: navier_stokes.cpp:62
MoFEM::SnesRhs
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:27
Range
SnesCtx
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
NavierStokesExample::setupAlgebraicStructures
MoFEMErrorCode setupAlgebraicStructures()
[Setup discrete manager]
Definition: navier_stokes.cpp:381
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:217
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
NavierStokesExample::numSteps
int numSteps
Definition: navier_stokes.cpp:53
NavierStokesExample::setupDiscreteManager
MoFEMErrorCode setupDiscreteManager()
[Define finite elements]
Definition: navier_stokes.cpp:357
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
NavierStokesExample::solveProblem
MoFEMErrorCode solveProblem()
[Setup SNES]
Definition: navier_stokes.cpp:556
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
NavierStokesExample::postProcess
MoFEMErrorCode postProcess()
[Solve problem]
Definition: navier_stokes.cpp:649
NavierStokesExample::fePostProcPtr
boost::shared_ptr< PostProcVol > fePostProcPtr
Definition: navier_stokes.cpp:79
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
NavierStokesExample::isHoGeometry
PetscBool isHoGeometry
Definition: navier_stokes.cpp:43
NavierStokesExample::feRhsPtr
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhsPtr
Definition: navier_stokes.cpp:71
NavierStokesExample::velocityScale
double velocityScale
Definition: navier_stokes.cpp:47
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
NavierStokesExample::numHoLevels
int numHoLevels
Definition: navier_stokes.cpp:41
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::check_finite_element
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
NavierStokesExample::lambdaStep
double lambdaStep
Definition: navier_stokes.cpp:54
MoFEM::SnesMat
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:139
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:590
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
NavierStokesExample::isStokesFlow
PetscBool isStokesFlow
Definition: navier_stokes.cpp:51
EshelbianPlasticity::U
@ U
Definition: EshelbianContact.cpp:193
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
NavierStokesExample::commonData
boost::shared_ptr< NavierStokesElement::CommonData > commonData
Definition: navier_stokes.cpp:64
NavierStokesElement::setNavierStokesOperators
static MoFEMErrorCode setNavierStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Navier-Stokes equations.
Definition: NavierStokesElement.cpp:22
NavierStokesElement::setStokesOperators
static MoFEMErrorCode setStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Stokes equations.
Definition: NavierStokesElement.cpp:79
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1127
MoFEM::DMMoFEMSNESSetFunction
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMoFEM.cpp:722
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
NavierStokesExample::setupSNES
MoFEMErrorCode setupSNES()
[Setup element instances]
Definition: navier_stokes.cpp:520