10using namespace boost::numeric;
13static char help[] =
"-order approximation order\n"
64 tD(
i,
j,
k,
l) *= coefficient;
106 if (row_side == col_side && row_type == col_type) {
139 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
141 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
142 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
143 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2));
152 double vol = getVolume();
155 auto t_w = getFTensor0IntegrationWeight();
163 const double a = t_w * vol;
166 for (
int rr = 0; rr !=
nbRows / 3; ++rr) {
169 auto t_m = get_tensor2(
K, 3 * rr, 0);
175 for (
int cc = 0; cc !=
nbCols / 3; ++cc) {
179 a * (
tD(
i,
j,
k,
l) * (t_row_diff_base(
j) * t_col_diff_base(
l)));
209 const int *row_indices = &*row_data.
getIndices().data().begin();
211 const int *col_indices = &*col_data.
getIndices().data().begin();
212 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
213 : getFEMethod()->snes_B;
216 &*
K.data().begin(), ADD_VALUES);
223 &*
K.data().begin(), ADD_VALUES);
254 nF.resize(nb_dofs,
false);
258 const int nb_gauss_pts = data.
getN().size1();
271 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
273 double w = 0.5 * t_w;
279 for (
int bb = 0; bb != nb_dofs / 3; ++bb) {
306 const Range &fix_second_node)
315 std::set<int> set_fix_dofs;
318 if (dit->get()->getDofCoeffIdx() == 2) {
320 set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
325 if (dit->get()->getDofCoeffIdx() == 1) {
326 set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
331 set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
335 std::vector<int> fix_dofs(set_fix_dofs.size());
337 std::copy(set_fix_dofs.begin(), set_fix_dofs.end(), fix_dofs.begin());
339 CHKERR MatAssemblyBegin(
ksp_B, MAT_FINAL_ASSEMBLY);
348 CHKERR MatZeroRowsColumns(
ksp_B, fix_dofs.size(), &fix_dofs[0], 1,
x,
390int main(
int argc,
char *argv[]) {
394 "-pc_factor_mat_solver_type mumps \n"
395 "-mat_mumps_icntl_20 0 \n"
399 if (!
static_cast<bool>(ifstream(
param_file))) {
400 std::ofstream file(
param_file.c_str(), std::ios::ate);
401 if (file.is_open()) {
410 moab::Core mb_instance;
411 moab::Interface &moab = mb_instance;
422 PetscBool flg_test = PETSC_FALSE;
423 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"SimpleElasticProblem",
431 CHKERR PetscOptionsBool(
"-test",
"if true is ctest",
"", flg_test,
432 &flg_test, PETSC_NULL);
434 ierr = PetscOptionsEnd();
442 Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
447 BRICK_PRESSURE_FACES = 3,
453 int id =
bit->getMeshsetId();
455 if (
id == FIX_BRICK_FACES) {
461 CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 0,
false, adj_ents,
462 moab::Interface::UNION);
464 CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 1,
false, adj_ents,
465 moab::Interface::UNION);
466 fix_faces.merge(adj_ents);
467 }
else if (
id == FIX_NODES) {
472 }
else if (
id == BRICK_PRESSURE_FACES) {
474 meshset, 2, pressure_faces,
true);
479 meshset, 0, fix_second_node,
true);
516 boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
518 boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
522 elastic_fe->getRuleHook =
VolRule();
523 pressure_fe->getRuleHook =
FaceRule();
527 elastic_fe->getOpPtrVector().push_back(
new OpK());
529 pressure_fe->getOpPtrVector().push_back(
new OpPressure());
531 boost::shared_ptr<FEMethod> fix_dofs_fe(
534 boost::shared_ptr<FEMethod> null_fe;
538 dm, simple_interface->
getDomainFEName(), elastic_fe, null_fe, null_fe);
553 CHKERR DMCreateGlobalVector(dm, &x);
557 CHKERR VecDuplicate(x, &f);
561 elastic_fe->ksp_B =
A;
562 fix_dofs_fe->ksp_B =
A;
566 fix_dofs_fe->ksp_f = f;
567 pressure_fe->ksp_f = f;
583 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
587 CHKERR KSPSetFromOptions(solver);
590 CHKERR KSPSetOperators(solver,
A,
A);
596 CHKERR KSPSolve(solver, f, x);
599 CHKERR KSPDestroy(&solver);
603 VecView(x, PETSC_VIEWER_STDOUT_WORLD);
610 auto get_post_proc_ele = [&]() {
611 auto jac_ptr = boost::make_shared<MatrixDouble>();
612 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
613 auto det_ptr = boost::make_shared<VectorDouble>();
615 auto post_proc_ele = boost::make_shared<
619 post_proc_ele->getOpPtrVector().push_back(
621 post_proc_ele->getOpPtrVector().push_back(
623 post_proc_ele->getOpPtrVector().push_back(
626 auto u_ptr = boost::make_shared<MatrixDouble>();
627 auto grad_ptr = boost::make_shared<MatrixDouble>();
628 post_proc_ele->getOpPtrVector().push_back(
630 post_proc_ele->getOpPtrVector().push_back(
635 post_proc_ele->getOpPtrVector().push_back(
639 post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),
644 {{
"GRAD", grad_ptr}},
650 return post_proc_ele;
653 if (
auto post_proc = get_post_proc_ele()) {
657 CHKERR post_proc->writeFile(
"out.h5m");
661 if (flg_test == PETSC_TRUE) {
662 const double x_vec_norm_const = 0.4;
667 CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
668 if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const >
671 "test failed (nrm 2 %6.4e)", norm_check);
const std::string default_options
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
ApplyDirichletBc(const Range &fix_faces, const Range &fix_nodes, const Range &fix_second_node)
MoFEMErrorCode postProcess()
function is run at the end of loop
Set integration rule to boundary elements.
int operator()(int, int, int p) 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 field gradients at integration pts for scalar filed rank 0, i.e. vector field.
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.
OpPressure(const double pressure_val=1)
int operator()(int, int, int p) const