24static char help[] =
"...\n\n";
78 "out_level_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
85 boost::shared_ptr<PostProcEle>
postProc;
168 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"ESSENTIAL",
171 auto &bc_map = bc_mng->getBcMapByBlockName();
172 boundaryMarker = boost::make_shared<std::vector<char unsigned>>();
173 for (
auto b : bc_map) {
174 if (std::regex_match(b.first, std::regex(
"(.*)ESSENTIAL(.*)"))) {
176 for (
int i = 0;
i != b.second->bcMarkers.size(); ++
i) {
177 (*boundaryMarker)[
i] |= b.second->bcMarkers[
i];
188 auto add_domain_base_ops = [&](
auto &pipeline) {
189 auto det_ptr = boost::make_shared<VectorDouble>();
190 auto jac_ptr = boost::make_shared<MatrixDouble>();
191 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
198 auto add_domain_lhs_ops = [&](
auto &pipeline) {
201 "U",
"U", [](
double,
double,
double) ->
double {
return 1; }));
211 auto add_domain_rhs_ops = [&](
auto &pipeline) {
213 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
214 auto dot2_u_at_gauss_pts = boost::make_shared<VectorDouble>();
216 "U", grad_u_at_gauss_pts));
220 "U", grad_u_at_gauss_pts,
221 [](
double,
double,
double) ->
double {
return 1; }));
223 "U", dot2_u_at_gauss_pts,
224 [](
const double,
const double,
const double) {
return wave_speed2; }));
228 auto add_boundary_base_ops = [&](
auto &pipeline) {};
230 auto add_lhs_base_ops = [&](
auto &pipeline) {
233 "U",
"U", [](
const double,
const double,
const double) {
return 1; }));
237 auto add_rhs_base_ops = [&](
auto &pipeline) {
239 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
240 auto boundary_function = [&](
const double x,
const double y,
244 const double t = fe_rhs->ts_t;
245 if ((
t <= 0.5) && (x < 0.) && (y > -1. / 3) && (y < 1. / 3))
246 return sin(4 * M_PI *
t);
253 [](
const double,
const double,
const double) {
return 1; }));
259 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
260 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
261 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
262 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
264 add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
265 add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
266 add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
267 add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
278 auto create_post_process_element = [&]() {
279 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
281 auto det_ptr = boost::make_shared<VectorDouble>();
282 auto jac_ptr = boost::make_shared<MatrixDouble>();
283 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
285 post_proc_fe->getOpPtrVector().push_back(
287 post_proc_fe->getOpPtrVector().push_back(
290 auto u_ptr = boost::make_shared<VectorDouble>();
291 post_proc_fe->getOpPtrVector().push_back(
296 post_proc_fe->getOpPtrVector().push_back(
298 new OpPPMap(post_proc_fe->getPostProcMesh(),
299 post_proc_fe->getMapGaussPts(),
312 auto set_time_monitor = [&](
auto dm,
auto solver) {
314 boost::shared_ptr<Monitor> monitor_ptr(
315 new Monitor(dm, create_post_process_element()));
316 boost::shared_ptr<ForcesAndSourcesCore> null;
318 monitor_ptr, null, null);
322 auto set_fieldsplit_preconditioner = [&](
auto solver) {
325 CHKERR TSGetSNES(solver, &snes);
327 CHKERR SNESGetKSP(snes, &ksp);
329 CHKERR KSPGetPC(ksp, &pc);
330 PetscBool is_pcfs = PETSC_FALSE;
331 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
333 if (is_pcfs == PETSC_TRUE) {
334 auto bc_mng = mField.getInterface<
BcManager>();
335 auto name_prb =
simple->getProblemName();
336 auto is_all_bc = bc_mng->getBlockIS(name_prb,
"ESSENTIAL",
"U", 0, 0);
338 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
340 <<
"Field split block size " << is_all_bc_size;
341 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
347 auto dm =
simple->getDM();
351 auto ts = pipeline_mng->createTSIM2();
355 CHKERR TS2SetSolution(ts,
D, DD);
356 CHKERR set_time_monitor(dm, ts);
357 CHKERR TSSetFromOptions(ts);
358 CHKERR set_fieldsplit_preconditioner(ts);
362 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
363 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
392int main(
int argc,
char *argv[]) {
399 auto core_log = logging::core::get();
408 DMType dm_name =
"DMMOFEM";
412 moab::Core mb_instance;
413 moab::Interface &moab = mb_instance;
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#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 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 DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
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.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_TT > OpCalculateScalarFieldValuesDotDot
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
constexpr double t
plate stiffness
static constexpr int approx_order
Simple interface for fast problem set-up.
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)
structure for User Loop Methods on finite elements
default operator for TRI element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
static constexpr int saveEveryNthStep
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode preProcess()
function is run at the beginning of loop
WaveEquation(MoFEM::Interface &m_field)
MoFEM::Interface & mField
MoFEMErrorCode readMesh()
MoFEMErrorCode setupProblem()
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode initialCondition()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode assembleSystem()
MoFEMErrorCode solveSystem()
MoFEMErrorCode outputResults()
MoFEMErrorCode runProgram()
MoFEMErrorCode setIntegrationRules()
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
constexpr double wave_speed2
constexpr int SPACE_DIM
[Define dimension]
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass