The example shows how to solve the linear elastic problem.
using namespace boost::numeric;
static char help[] =
"-order approximation order\n"
"\n";
symm) {
tD(
i,
j,
k,
l) *= coefficient;
}
if (row_side == col_side && row_type == col_type) {
} else {
}
}
protected:
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
&
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
&
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2));
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
const double a = t_w * vol;
for (
int rr = 0; rr !=
nbRows / 3; ++rr) {
auto t_m = get_tensor2(
K, 3 * rr, 0);
for (
int cc = 0; cc !=
nbCols / 3; ++cc) {
a * (
tD(
i,
j,
k,
l) * (t_row_diff_base(
j) * t_col_diff_base(
l)));
++t_col_diff_base;
++t_m;
}
++t_row_diff_base;
}
++t_w;
}
}
const int *row_indices = &*row_data.
getIndices().data().begin();
const int *col_indices = &*col_data.
getIndices().data().begin();
Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
: getFEMethod()->snes_B;
&*
K.data().begin(), ADD_VALUES);
&*
K.data().begin(), ADD_VALUES);
}
}
};
if (nb_dofs == 0)
nF.resize(nb_dofs,
false);
const int nb_gauss_pts = data.
getN().size1();
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
for (int bb = 0; bb != nb_dofs / 3; ++bb) {
++t_nf;
++t_base;
}
++t_w;
}
}
};
const Range &fix_second_node)
}
std::set<int> set_fix_dofs;
if (dit->get()->getDofCoeffIdx() == 2) {
set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
}
}
if (dit->get()->getDofCoeffIdx() == 1) {
set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
}
}
set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
}
}
std::vector<int> fix_dofs(set_fix_dofs.size());
std::copy(set_fix_dofs.begin(), set_fix_dofs.end(), fix_dofs.begin());
CHKERR MatZeroRowsColumns(
ksp_B, fix_dofs.size(), &fix_dofs[0], 1,
x,
}
};
int operator()(
int,
int,
int p)
const {
return 2 * (
p - 1); }
};
};
int main(
int argc,
char *argv[]) {
"-pc_type lu \n"
"-pc_factor_mat_solver_type mumps \n"
"-mat_mumps_icntl_20 0 \n"
"-ksp_monitor\n";
std::ofstream file(
param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file.close();
}
}
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
try {
PetscBool flg_test = PETSC_FALSE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"SimpleElasticProblem",
"none");
PETSC_NULL);
CHKERR PetscOptionsBool(
"-test",
"if true is ctest",
"", flg_test,
&flg_test, PETSC_NULL);
ierr = PetscOptionsEnd();
Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
enum MyBcTypes {
FIX_BRICK_FACES = 1,
FIX_NODES = 2,
BRICK_PRESSURE_FACES = 3,
FIX_NODES_Y_DIR = 4
};
int id =
bit->getMeshsetId();
if (id == FIX_BRICK_FACES) {
fix_faces, true);
CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 0,
false, adj_ents,
moab::Interface::UNION);
CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 1,
false, adj_ents,
moab::Interface::UNION);
fix_faces.merge(adj_ents);
} else if (id == FIX_NODES) {
fix_nodes, true);
} else if (id == BRICK_PRESSURE_FACES) {
meshset, 2, pressure_faces, true);
} else if (id ==
FIX_NODES_Y_DIR) {
meshset, 0, fix_second_node, true);
} else {
}
}
3);
DM dm;
"PRESSURE");
boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
elastic_fe->getRuleHook =
VolRule();
elastic_fe->getOpPtrVector().push_back(
new OpK());
pressure_fe->getOpPtrVector().push_back(
new OpPressure());
boost::shared_ptr<FEMethod> fix_dofs_fe(
boost::shared_ptr<FEMethod> null_fe;
null_fe);
CHKERR DMCreateGlobalVector(dm, &x);
elastic_fe);
KSP solver;
CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSolve(solver, f, x);
VecView(x, PETSC_VIEWER_STDOUT_WORLD);
auto get_post_proc_ele = [&]() {
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto post_proc_ele = boost::make_shared<
post_proc_ele->getOpPtrVector().push_back(
post_proc_ele->getOpPtrVector().push_back(
post_proc_ele->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
post_proc_ele->getOpPtrVector().push_back(
post_proc_ele->getOpPtrVector().push_back(
grad_ptr));
post_proc_ele->getOpPtrVector().push_back(
post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),
{},
{{"U", u_ptr}},
{{"GRAD", grad_ptr}},
{})
);
return post_proc_ele;
};
if (auto post_proc = get_post_proc_ele()) {
post_proc);
CHKERR post_proc->writeFile(
"out.h5m");
}
{
if (flg_test == PETSC_TRUE) {
const double x_vec_norm_const = 0.4;
double norm_check;
CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const >
1.e-10) {
"test failed (nrm 2 %6.4e)", norm_check);
}
}
}
}
return 0;
}
const std::string default_options
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_NUMEREDDOF_ROW_FOR_LOOP_(PROBLEMPTR, IT)
use with loops to iterate row DOFs
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
MoFEMErrorCode postProcess()
function is run at the end of loop
Set integration rule to boundary elements.
int operator()(int, int, int) const
const Problem * problemPtr
raw pointer to problem
virtual moab::Interface & get_moab()=0
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure for User Loop Methods on finite elements
default operator for TRI element
auto getFTensor1Normal()
get normal as tensor
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Simple interface for fast problem set-up.
MoFEMErrorCode buildProblem()
Build problem.
MoFEMErrorCode defineFiniteElements()
Define finite elements.
MoFEMErrorCode buildFiniteElements()
Build finite elements.
MoFEMErrorCode addDomainField(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_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode buildFields()
Build fields.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
const std::string getDomainFEName() const
Get the Domain FE Name.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate B^T D B operator.
bool isDiag
true if this block is on diagonal
FTensor::Ddg< double, 3, 3 > tD
int nbIntegrationPts
number of integration points
int nbRows
number of dofs on rows
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations for give operator.
int nbCols
number if dof on column
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble local entity block matrix.
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
int operator()(int, int, int) const