#include <stdlib.h>
#include <BasicFiniteElements.hpp>
static char help[] =
"...\n\n";
public:
boost::shared_ptr<MatrixDouble> field_grad_mat)
fieldGradMat(field_grad_mat) {}
const double area = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_field_grad = getFTensor1FromMat<2>(*fieldGradMat);
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
const double an = 1. / std::sqrt(1 + t_field_grad(
i) * t_field_grad(
i));
locLhs(rr, cc) += (t_row_diff_base(
i) * t_col_diff_base(
i)) * an *
a -
(t_field_grad(
i) * t_col_diff_base(
i)) *
t_field_grad(
i) * an * an * an *
a;
++t_col_diff_base;
}
++t_row_diff_base;
}
++t_w;
++t_field_grad;
}
}
private:
boost::shared_ptr<MatrixDouble> fieldGradMat;
};
public:
boost::shared_ptr<MatrixDouble> field_grad_mat)
fieldGradMat(field_grad_mat) {}
const double area = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_field_grad = getFTensor1FromMat<2>(*fieldGradMat);
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
const double an = 1. / std::sqrt(1 + t_field_grad(
i) * t_field_grad(
i));
nf[rr] += (t_diff_base(
i) * t_field_grad(
i)) * an *
a;
++t_base;
++t_diff_base;
}
++t_w;
++t_field_grad;
}
}
private:
boost::shared_ptr<MatrixDouble> fieldGradMat;
};
public:
private:
static double boundaryFunction(const double x, const double y,
const double z) {
return sin(2 * M_PI * (x + y));
}
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
Range boundaryEnts;
};
: mField(m_field) {}
}
}
}
};
}
auto get_ents_on_mesh_skin = [&]() {
Range body_ents;
Range skin_ents;
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
Range boundary_ents;
ParallelComm *pcomm =
pcomm->filter_pstatus(skin_ents, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
Range skin_verts;
boundary_ents.merge(skin_verts);
return boundary_ents;
};
auto mark_boundary_dofs = [&](Range &&skin_edges) {
auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
problem_manager->markDofs(
simple->getProblemName(),
ROW,
ProblemsManager::OR, skin_edges, *marker_ptr);
return marker_ptr;
};
}
auto add_domain_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
pipeline.push_back(new OpSetHOWeightsOnFace());
};
auto add_domain_lhs_ops = [&](auto &pipeline) {
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
pipeline.push_back(
new OpCalculateScalarFieldGradient<2>("U", grad_u_at_gauss_pts));
pipeline.push_back(
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_rhs_ops = [&](auto &pipeline) {
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
pipeline.push_back(
new OpCalculateScalarFieldGradient<2>("U", grad_u_at_gauss_pts));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_boundary_base_ops = [&](auto &pipeline) {};
auto add_lhs_base_ops = [&](auto &pipeline) {
"U", "U", [](const double, const double, const double) { return 1; }));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_rhs_base_ops = [&](auto &pipeline) {
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
"U", u_at_gauss_pts,
[](const double, const double, const double) { return 1; }));
pipeline.push_back(new OpUnSetBc("U"));
};
add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
}
auto set_fieldsplit_preconditioner = [&](auto snes) {
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto name_prb =
simple->getProblemName();
SmartPetscObj<IS> is_all_bc;
int is_all_bc_size;
CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
<< "Field split block size " << is_all_bc_size;
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc);
}
};
SmartPetscObj<Vec> global_rhs, global_solution;
auto solver = pipeline_mng->createSNES();
CHKERR SNESSetFromOptions(solver);
CHKERR set_fieldsplit_preconditioner(solver);
CHKERR SNESSolve(solver, global_rhs, global_solution);
SCATTER_REVERSE);
}
auto post_proc = boost::make_shared<PostProcEle>(
mField);
CHKERR post_proc->generateReferenceElementMesh();
CHKERR post_proc->addFieldValuesPostProc(
"U");
CHKERR post_proc->writeFile(
"out_result.h5m");
}
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
CHKERR minimal_surface_problem.runProgram();
}
return 0;
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
int main(int argc, char *argv[])
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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 DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBoundaryTimeScalarField
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode outputResults()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode runProgram()
MoFEMErrorCode solveSystem()
MoFEM::Interface & mField
static double boundaryFunction(const double x, const double y, const double z)
MinimalSurfaceEqn(MoFEM::Interface &m_field)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode setupProblem()
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.
default operator for TRI element
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Integrate the domain residual vector (RHS)
Integrate the domain tangent matrix (LHS)