#include <BasicFiniteElements.hpp>
using namespace boost::numeric;
using namespace std;
static char help[] =
"Navier-Stokes Example\n";
commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
}
}
private:
PetscBool isPartitioned = PETSC_FALSE;
int orderVelocity = 2;
int orderPressure = 1;
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;
int step;
Range solidFaces;
DM dM;
SNES snes;
boost::shared_ptr<NavierStokesElement::CommonData> commonData;
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<PostProcVolumeOnRefinedMesh> fePostProcPtr;
boost::shared_ptr<PostProcFaceOnRefinedMesh> fePostProcDragPtr;
};
CHKERR defineFiniteElements();
CHKERR setupDiscreteManager();
CHKERR setupAlgebraicStructures();
CHKERR setupElementInstances();
}
PetscBool flg_mesh_file;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"NAVIER_STOKES problem",
"none");
CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
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",
"",
isDiscontPressure, &isDiscontPressure, PETSC_NULL);
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)");
}
if (isPartitioned == PETSC_TRUE) {
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
} else {
const char *option;
option = "";
}
CHKERR mField.rebuild_database();
bitLevel.set(0);
CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(0, 3,
bitLevel);
}
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) {
"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);
CHKERR mField.get_moab().get_adjacencies(
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) {
"Cannot find a fluid block adjacent to a given solid face");
}
}
}
}
if (isDiscontPressure) {
} else {
}
3);
CHKERR mField.add_ents_to_field_by_dim(0, 3,
"VELOCITY");
CHKERR mField.add_ents_to_field_by_dim(0, 3,
"PRESSURE");
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);
if (!isDiscontPressure) {
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--) {
if (isPartitioned)
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
levels[ll]);
CHKERR mField.set_field_order(levels[ll],
"VELOCITY",
orderVelocity + add_order);
if (!isDiscontPressure)
CHKERR mField.set_field_order(levels[ll],
"PRESSURE",
orderPressure + add_order);
else
CHKERR mField.set_field_order(levels[ll].subset_by_type(MBTET),
"PRESSURE", orderPressure + add_order);
++add_order;
}
}
CHKERR mField.set_field_order(0, MBVERTEX,
"MESH_NODE_POSITIONS", 1);
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);
if (isPartitioned)
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(edges);
CHKERR mField.set_field_order(edges,
"MESH_NODE_POSITIONS", 2);
}
if (isPartitioned) {
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);
}
"PRESSURE", "MESH_NODE_POSITIONS");
"MESH_NODE_POSITIONS", 2, &solidFaces);
CHKERR mField.build_finite_elements();
CHKERR mField.build_adjacencies(bitLevel);
}
DMType dm_name = "DM_NAVIER_STOKES";
CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
CHKERR DMSetType(dM, dm_name);
if (mField.check_finite_element("PRESSURE_FE"))
}
CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
}
feLhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
feRhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
feDragPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
feDragSidePtr = boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>(
true);
if (isStokesFlow) {
feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
} else {
feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
}
"NAVIER_STOKES", "VELOCITY",
"PRESSURE", commonData);
dirichletBcPtr = boost::shared_ptr<DirichletDisplacementBc>(
dirichletBcPtr->snes_x =
D;
NULL, "VELOCITY");
fePostProcPtr = boost::make_shared<PostProcVolumeOnRefinedMesh>(mField);
CHKERR fePostProcPtr->generateReferenceElementMesh();
true);
CHKERR fePostProcPtr->addFieldValuesPostProc(
"VELOCITY");
CHKERR fePostProcPtr->addFieldValuesPostProc(
"PRESSURE");
CHKERR fePostProcPtr->addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
CHKERR fePostProcPtr->addFieldValuesGradientPostProc(
"VELOCITY");
fePostProcDragPtr = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
CHKERR fePostProcDragPtr->generateReferenceElementMesh();
CHKERR fePostProcDragPtr->addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
fePostProcDragPtr, feDragSidePtr, "NAVIER_STOKES", "VELOCITY", "PRESSURE",
commonData);
}
boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
neumannForces.begin();
for (; mit != neumannForces.end(); mit++) {
&mit->second->getLoopFe(), NULL, NULL);
}
boost::shared_ptr<FEMethod> null_fe;
null_fe);
dirichletBcPtr);
null_fe);
null_fe);
dirichletBcPtr);
CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
CHKERR SNESSetFromOptions(snes);
}
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");
};
pressureScale = viscosity * velocityScale / lengthScale;
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",
velocityScale);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Pressure scale: %6.4e\n",
pressureScale);
if (isStokesFlow) {
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Re number: 0 (Stokes flow)\n");
} else {
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Re number: %6.4e\n",
density / viscosity * velocityScale * lengthScale);
}
lambdaStep = 1.0 / numSteps;
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 {
reNumber = density / viscosity * velocityScale * lengthScale *
lambda;
for (
auto &
bit : commonData->setOfBlocksData) {
bit.second.inertiaCoef = reNumber;
bit.second.viscousCoef = 1.0;
}
}
PETSC_COMM_WORLD,
"Step: %d | Lambda: %6.4e | Inc: %6.4e | Re number: %6.4e \n", step,
lambda, lambdaStep, reNumber);
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 scale_problem(velocityScale, lengthScale, pressureScale);
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->postProcMesh.write_file(out_file_name.c_str(),
"MOAB",
"PARALLEL=WRITE_PART");
}
CHKERR VecZeroEntries(commonData->pressureDragForceVec);
CHKERR VecZeroEntries(commonData->shearDragForceVec);
CHKERR VecZeroEntries(commonData->totalDragForceVec);
auto get_vec_data = [&](auto vec, std::array<double, 3> &data) {
const double *array;
CHKERR VecGetArrayRead(vec, &array);
if (mField.get_comm_rank() == 0) {
}
CHKERR VecRestoreArrayRead(vec, &array);
};
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) {
"Pressure drag force: [ %6.4e, %6.4e, %6.4e ]",
pressure_drag[0], pressure_drag[1], pressure_drag[2]);
"Shear drag force: [ %6.4e, %6.4e, %6.4e ]", shear_drag[0],
shear_drag[1], shear_drag[2]);
"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());
CHKERR fePostProcDragPtr->postProcMesh.write_file(
out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
}
}
int main(
int argc,
char *argv[]) {
try {
}
return 0;
}
#define MOFEM_LOG_C(channel, severity, format,...)
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto smartCreateDMVector
Get smart vector from DM.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
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.
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
auto smartCreateDMMatrix
Get smart matrix from DM.
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreOnSideSwitch< 0 > VolumeElementForcesAndSourcesCoreOnSide
Volume element used to integrate on skeleton.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#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.
FTensor::Index< 'i', SPACE_DIM > i
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
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.
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.
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
const double D
diffusivity
int main(int argc, char *argv[])
[Post process]
Set Dirichlet boundary conditions on displacements.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Set integration rule to volume elements.
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.
static MoFEMErrorCode setPostProcDragOperators(boost::shared_ptr< PostProcFaceOnRefinedMesh > 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.
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.
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.
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.
MoFEMErrorCode setupElementInstances()
[Setup algebraic structures]
MoFEMErrorCode setupDiscreteManager()
[Define finite elements]
MoFEMErrorCode defineFiniteElements()
[Setup fields]
MoFEMErrorCode setupFields()
[Find blocks]
MoFEMErrorCode setupSNES()
[Setup element instances]
MoFEMErrorCode runProblem()
[Example Navier Stokes]
MoFEMErrorCode setupAlgebraicStructures()
[Setup discrete manager]
MoFEMErrorCode findBlocks()
[Read input]
MoFEMErrorCode solveProblem()
[Setup SNES]
MoFEMErrorCode postProcess()
[Solve problem]
MoFEMErrorCode readInput()
[Run problem]