\[ \begin{aligned} \frac{\partial u(\mathbf{x}, t)}{\partial t}-\Delta u(\mathbf{x}, t) &=f(\mathbf{x}, t) & & \forall \mathbf{x} \in \Omega, t \in(0, T), \\ u(\mathbf{x}, 0) &=u_{0}(\mathbf{x}) & & \forall \mathbf{x} \in \Omega, \\ u(\mathbf{x}, t) &=g(\mathbf{x}, t) & & \forall \mathbf{x} \in \partial \Omega, t \in(0, T). \end{aligned} \]
#include <stdlib.h>
#include <cmath>
#include <BasicFiniteElements.hpp>
static char help[] =
"...\n\n";
public:
private:
static double sourceTermFunction(const double x, const double y,
const double z, const double t) {
return 0.1 * pow(M_E, -M_PI * M_PI * t) * sin(1. * M_PI * x) *
sin(2. * M_PI * y);
}
static double boundaryFunction(const double x, const double y, const double z,
const double t) {
return abs(0.1 * pow(M_E, -M_PI * M_PI * t) * sin(2. * M_PI * x) *
sin(3. * M_PI * y));
}
const double initU = 0.1;
Simple *simpleInterface;
MPI_Comm mpiComm;
const int mpiRank;
SmartPetscObj<DM> dM;
SmartPetscObj<TS> tsSolver;
std::string domainField;
int oRder;
boost::shared_ptr<std::vector<bool>> boundaryMarker;
boost::shared_ptr<FaceEle> stiffTangentLhsMatrixPipeline;
boost::shared_ptr<FaceEle> stiffFunctionRhsVectorPipeline;
boost::shared_ptr<EdgeEle> boundaryLhsMatrixPipeline;
boost::shared_ptr<EdgeEle> boundaryRhsVectorPipeline;
boost::shared_ptr<DataAtGaussPoints> previousUpdate;
boost::shared_ptr<VectorDouble> fieldValuePtr;
boost::shared_ptr<MatrixDouble> fieldGradPtr;
boost::shared_ptr<VectorDouble> fieldDotPtr;
boost::shared_ptr<PostProcFaceOnRefinedMesh> postProcFace;
Range boundaryEntitiesForFieldsplit;
};
: domainField("U"), mField(m_field), mpiComm(mField.get_comm()),
mpiRank(mField.get_comm_rank()) {
stiffTangentLhsMatrixPipeline =
boost::shared_ptr<FaceEle>(
new FaceEle(mField));
stiffFunctionRhsVectorPipeline =
boost::shared_ptr<FaceEle>(
new FaceEle(mField));
boundaryLhsMatrixPipeline = boost::shared_ptr<EdgeEle>(
new EdgeEle(mField));
boundaryRhsVectorPipeline = boost::shared_ptr<EdgeEle>(
new EdgeEle(mField));
previousUpdate =
boost::shared_ptr<DataAtGaussPoints>(new DataAtGaussPoints());
fieldValuePtr = boost::shared_ptr<VectorDouble>(previousUpdate,
&previousUpdate->fieldValue);
fieldGradPtr = boost::shared_ptr<MatrixDouble>(previousUpdate,
&previousUpdate->fieldGrad);
fieldDotPtr = boost::shared_ptr<VectorDouble>(previousUpdate,
&previousUpdate->fieldDot);
}
}
}
auto stiff_tangent_rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
auto stiff_function_rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
auto boundary_lhs_rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
auto boundary_rhs_rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
}
Range inner_surface;
if (!inner_surface.empty()) {
Range inner_surface_verts;
inner_surface_verts, false);
}
}
}
auto get_ents_on_mesh_skin = [&]() {
Range boundary_entities;
std::string entity_name = it->getName();
if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
boundary_entities, true);
}
}
Range boundary_vertices;
boundary_vertices, true);
boundary_entities.merge(boundary_vertices);
return boundary_entities;
};
auto mark_boundary_dofs = [&](Range &&skin_edges) {
auto marker_ptr = boost::make_shared<std::vector<bool>>();
skin_edges, *marker_ptr);
return marker_ptr;
};
}
{
new OpSetInvJacH1ForFace(
invJac));
}
{
new OpSetInvJacH1ForFace(
invJac));
}
{
}
}
boost::shared_ptr<ForcesAndSourcesCore> null;
{
postProcFace = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
monitor_ptr, null, null);
}
SmartPetscObj<Vec> global_solution;
SCATTER_FORWARD);
}
}
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
CHKERR heat_problem.runProgram();
}
return 0;
}
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCoreSwitch< EdgeElementForcesAndSourcesCore::NO_HO_GEOMETRY > EdgeEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#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.
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
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 DMMoFEMTSSetIJacobian(DM dm, 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 TS Jacobian evaluation function
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
#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.
int main(int argc, char *argv[])
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
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.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
DeprecatedCoreInterface Interface
boost::shared_ptr< FaceEle > stiffFunctionRhsVectorPipeline
boost::shared_ptr< EdgeEle > boundaryRhsVectorPipeline
SmartPetscObj< TS > tsSolver
boost::shared_ptr< std::vector< bool > > boundaryMarker
MoFEM::Interface & mField
MoFEMErrorCode solveSystem()
MoFEMErrorCode outputResults()
boost::shared_ptr< FaceEle > stiffTangentLhsMatrixPipeline
boost::shared_ptr< EdgeEle > boundaryLhsMatrixPipeline
Range boundaryEntitiesForFieldsplit
boost::shared_ptr< VectorDouble > fieldDotPtr
boost::shared_ptr< VectorDouble > fieldValuePtr
HeatEquation(MoFEM::Interface &m_field)
boost::shared_ptr< DataAtGaussPoints > previousUpdate
MoFEMErrorCode setupProblem()
static double sourceTermFunction(const double x, const double y, const double z, const double t)
MoFEMErrorCode setIntegrationRules()
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcFace
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode runProgram()
MoFEMErrorCode initialCondition()
static double boundaryFunction(const double x, const double y, const double z, const double t)
boost::shared_ptr< MatrixDouble > fieldGradPtr
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =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.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
[Push operators to pipeline]