v0.14.0
Loading...
Searching...
No Matches
navier_stokes.cpp

Example of viscous fluid flow problem

/** \file navier_stokes.cpp
* \example navier_stokes.cpp
*
* Example of viscous fluid flow problem
*
*/
using PostProcVol = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
using PostProcFace = PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore>;
using namespace boost::numeric;
using namespace MoFEM;
using namespace std;
static char help[] = "Navier-Stokes Example\n";
//! [Example Navier Stokes]
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;
BitRefLevel bitLevel;
DM dM;
SNES snes;
boost::shared_ptr<NavierStokesElement::CommonData> commonData;
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;
};
//! [Example Navier Stokes]
//! [Run problem]
}
//! [Run problem]
//! [Read input]
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;
option = "PARALLEL=READ_PART;"
"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,
}
//! [Read input]
//! [Find blocks]
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);
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) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Cannot find a fluid block adjacent to a given solid face");
}
}
}
}
//! [Find blocks]
//! [Setup fields]
} 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",
orderVelocity + add_order);
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);
// 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");
}
"MESH_NODE_POSITIONS");
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
}
//! [Setup fields]
//! [Define finite elements]
// 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);
// setup elements for loading
// build finite elements
// build adjacencies between elements and degrees of freedom
}
//! [Define finite elements]
//! [Setup discrete manager]
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)
CHKERR DMMoFEMAddElement(dM, "NAVIER_STOKES");
if (mField.check_finite_element("PRESSURE_FE"))
CHKERR DMMoFEMAddElement(dM, "PRESSURE_FE");
// setup the DM
CHKERR DMSetUp(dM);
}
//! [Setup discrete manager]
//! [Setup algebraic structures]
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);
}
//! [Setup algebraic structures]
//! [Setup element instances]
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->methodsOp.push_back(new NavierStokesElement::LoadScale());
// 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(
new OpCalculateVectorFieldGradient<3, 3>("VELOCITY", grad_ptr));
fePostProcPtr->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PRESSURE", p_ptr));
fePostProcPtr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
fePostProcPtr->getOpPtrVector().push_back(
new OpPPMap(
fePostProcPtr->getPostProcMesh(), fePostProcPtr->getMapGaussPts(),
{{"PRESSURE", p_ptr}},
{{"VELOCITY", v_ptr}, {"MESH_NODE_POSITIONS", pos_ptr}},
{{"VELOCITY_GRAD", grad_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);
}
//! [Setup element instances]
//! [Setup SNES]
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);
}
//! [Setup SNES]
//! [Solve problem]
auto scale_problem = [&](double U, double L, double P) {
"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++;
}
}
//! [Solve problem]
//! [Post process]
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;
CHKERR VecGetArrayRead(vec, &array);
if (mField.get_comm_rank() == 0) {
for (int i : {0, 1, 2})
data[i] = array[i];
}
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) {
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());
}
}
//! [Post process]
//! [Main function]
int main(int argc, char *argv[]) {
const char param_file[] = "param_file.petsc";
// Initialise MoFEM
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;
}
//! [Main function]
#define DM_NO_ELEMENT
Definition DMMoFEM.hpp:10
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
static char help[]
int main()
PostProcEleByDim< SPACE_DIM >::PostProcFace PostProcFace
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ L2
field with C-1 continuity
Definition definitions.h:88
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ BLOCKSET
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:483
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:704
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
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:745
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
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:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1003
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:988
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition DMMoFEM.cpp:1080
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition DMMoFEM.cpp:521
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:532
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
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.
#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.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
char out_file_name[255]
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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
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
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore > PostProcVol
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Set Dirichlet boundary conditions on displacements.
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.
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.
Managing BitRefLevels.
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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
Add field.
virtual int get_comm_rank() const =0
virtual MoFEMErrorCode rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
Core (interface) class.
Definition Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
Deprecated interface functions.
Basic algebra on fields.
Definition FieldBlas.hpp:21
Get value at integration points for scalar field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
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 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 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.
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.
[Example Navier Stokes]
MoFEMErrorCode setupElementInstances()
[Setup algebraic structures]
MoFEMErrorCode setupDiscreteManager()
[Define finite elements]
boost::shared_ptr< FaceElementForcesAndSourcesCore > feDragPtr
MoFEMErrorCode defineFiniteElements()
[Setup fields]
boost::shared_ptr< PostProcFace > fePostProcDragPtr
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhsPtr
MoFEMErrorCode setupFields()
[Find blocks]
MoFEMErrorCode setupSNES()
[Setup element instances]
MoFEMErrorCode runProblem()
[Example Navier Stokes]
SmartPetscObj< Vec > F
boost::shared_ptr< NavierStokesElement::CommonData > commonData
boost::shared_ptr< DirichletDisplacementBc > dirichletBcPtr
MoFEMErrorCode setupAlgebraicStructures()
[Setup discrete manager]
NavierStokesExample(MoFEM::Interface &m_field)
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > feDragSidePtr
MoFEMErrorCode findBlocks()
[Read input]
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhsPtr
boost::ptr_map< std::string, NeumannForcesSurface > neumannForces
MoFEM::Interface & mField
boost::shared_ptr< PostProcVol > fePostProcPtr
MoFEMErrorCode solveProblem()
[Setup SNES]
MoFEMErrorCode postProcess()
[Solve problem]
MoFEMErrorCode readInput()
[Run problem]
SmartPetscObj< Vec > D
SmartPetscObj< Mat > Aij