using PostProcVol = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
using PostProcFace = PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore>;
using namespace boost::numeric;
static char help[] =
"Navier-Stokes Example\n";
commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
}
}
private:
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;
};
}
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?",
"",
CHKERR PetscOptionsInt(
"-my_order_u",
"approximation orderVelocity",
"",
CHKERR PetscOptionsInt(
"-my_order_p",
"approximation orderPressure",
"",
CHKERR PetscOptionsInt(
"-my_num_ho_levels",
"number of higher order boundary levels", "",
CHKERR PetscOptionsBool(
"-my_discont_pressure",
"discontinuous pressure",
"",
CHKERR PetscOptionsBool(
"-my_ho_geometry",
"use second order geometry",
"",
CHKERR PetscOptionsScalar(
"-my_velocity_scale",
"velocity scale",
"",
CHKERR PetscOptionsInt(
"-my_step_num",
"number of steps",
"",
numSteps,
ierr = PetscOptionsEnd();
if (flg_mesh_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
} else {
const char *option;
option = "";
}
}
if (
bit->getName().compare(0, 9,
"MAT_FLUID") == 0) {
const int id =
bit->getMeshsetId();
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].fluidViscosity = attributes[0];
commonData->setOfBlocksData[id].fluidDensity = attributes[1];
}
}
if (
bit->getName().compare(0, 9,
"INT_SOLID") == 0) {
const int id =
bit->getMeshsetId();
bit->getMeshset(), MBTRI,
commonData->setOfFacesData[
id].eNts,
true);
commonData->setOfFacesData[
id].eNts, 3,
true, tets,
moab::Interface::UNION);
tet =
Range(tets.front(), tets.front());
if (
bit.second.eNts.contains(tet)) {
bit.second.fluidViscosity;
commonData->setOfFacesData[id].fluidDensity =
bit.second.fluidDensity;
break;
}
}
if (
commonData->setOfFacesData[
id].fluidViscosity < 0) {
"Cannot find a fluid block adjacent to a given solid face");
}
}
}
}
} else {
}
3);
}
for (auto d : {1, 2, 3}) {
moab::Interface::UNION);
}
levels[ll] = subtract(ents, ents.subset_by_type(MBVERTEX));
}
levels[ll] = subtract(levels[ll], levels[ll - 1]);
}
int add_order = 1;
levels[ll]);
else
++add_order;
}
}
moab::Interface::UNION);
}
"VELOCITY");
"PRESSURE");
"MESH_NODE_POSITIONS");
}
"MESH_NODE_POSITIONS");
}
"PRESSURE", "MESH_NODE_POSITIONS");
}
DMType dm_name = "DM_NAVIER_STOKES";
}
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>(
true);
true);
feDragPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
feDragSidePtr = boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>(
true);
} else {
}
"NAVIER_STOKES", "VELOCITY",
NULL, "VELOCITY");
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>();
{{"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(
fePostProcDragPtr->getOpPtrVector().push_back(
fePostProcDragPtr->getPostProcMesh(),
fePostProcDragPtr->getMapGaussPts(),
{},
{{"MESH_NODE_POSITIONS", pos_ptr}},
{},
{}
)
);
fePostProcDragPtr, feDragSidePtr, "NAVIER_STOKES", "VELOCITY", "PRESSURE",
commonData);
}
boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
&mit->second->getLoopFe(), NULL, NULL);
}
boost::shared_ptr<FEMethod> null_fe;
null_fe);
null_fe);
null_fe);
}
auto scale_problem = [&](
double U,
double L,
double P) {
"MESH_NODE_POSITIONS");
mField,
"MESH_NODE_POSITIONS",
true,
true);
ent_method_on_10nodeTet.setNodes = false;
};
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Density: %6.4e\n",
density);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Velocity scale: %6.4e\n",
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Pressure scale: %6.4e\n",
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) {
bit.second.inertiaCoef = 0.0;
bit.second.viscousCoef = 1.0;
}
} else {
bit.second.viscousCoef = 1.0;
}
}
PETSC_COMM_WORLD,
"Step: %d | Lambda: %6.4e | Inc: %6.4e | Re number: %6.4e \n",
step,
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
}
{
std::ostringstream stm;
stm <<
"out_" <<
step <<
".h5m";
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
}
auto get_vec_data = [&](auto vec, std::array<double, 3> &data) {
const double *array;
CHKERR VecGetArrayRead(vec, &array);
}
CHKERR VecRestoreArrayRead(vec, &array);
};
std::array<double, 3> pressure_drag;
std::array<double, 3> shear_drag;
std::array<double, 3> total_drag;
"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";
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"out file %s\n",
}
}
int main(
int argc,
char *argv[]) {
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
return 0;
}
#define MOFEM_LOG_C(channel, severity, format,...)
static PetscErrorCode ierr
#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.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
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 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.
auto createDMVector(DM dm)
Get smart vector from DM.
auto createDMMatrix(DM dm)
Get smart matrix from 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)
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.
FTensor::Index< 'i', SPACE_DIM > i
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
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
PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore > PostProcFace
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Set Dirichlet boundary conditions on displacements.
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.
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.
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.
Volume finite element base.
Base volume element used to integrate on skeleton.
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.
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]
boost::shared_ptr< NavierStokesElement::CommonData > commonData
boost::shared_ptr< DirichletDisplacementBc > dirichletBcPtr
PetscBool isDiscontPressure
MoFEMErrorCode setupAlgebraicStructures()
[Setup discrete manager]
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]