#ifndef __NONLINEARPOISSON2D_HPP__
#define __NONLINEARPOISSON2D_HPP__
#include <stdlib.h>
#include <BasicFiniteElements.hpp>
typedef boost::function<double(const double, const double, const double)>
struct DataAtGaussPoints {
};
public:
std::string row_field_name, std::string col_field_name,
boost::shared_ptr<DataAtGaussPoints> &common_data,
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
}
EntityType col_type,
EntData &row_data,
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
locLhs.resize(nb_row_dofs, nb_col_dofs,
false);
const int nb_integration_points =
getGaussPts().size2();
auto t_field_grad = getFTensor1FromMat<2>(
commonData->fieldGrad);
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; cc++) {
locLhs(rr, cc) += (((1 + t_field * t_field) * t_row_diff_base(
i) *
(2.0 * t_field * t_field_grad(
i) *
t_row_diff_base(
i) * t_col_base)) *
++t_col_base;
++t_col_diff_base;
}
++t_row_base;
++t_row_diff_base;
}
++t_w;
++t_field;
++t_field_grad;
}
}
}
}
ADD_VALUES);
}
}
private:
};
public:
std::string field_name,
ScalarFunc source_term_function,
boost::shared_ptr<DataAtGaussPoints> &common_data,
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
if (nb_dofs) {
locRhs.resize(nb_dofs,
false);
const int nb_integration_points =
getGaussPts().size2();
auto t_field_grad = getFTensor1FromMat<2>(
commonData->fieldGrad);
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_dofs; rr++) {
t_grad_base(
i) * t_field_grad(
i) * (1 + t_field * t_field)) *
++t_base;
++t_grad_base;
}
++t_w;
++t_coords;
++t_field;
++t_field_grad;
}
}
}
}
CHKERR VecSetOption(
getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
}
}
private:
};
struct OpBoundaryTangentMatrix :
public OpEdgeEle {
public:
std::string col_field_name)
sYmm = true;
}
EntityType col_type,
EntData &row_data,
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
const double edge = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; cc++) {
++t_col_base;
}
++t_row_base;
}
++t_w;
}
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
}
}
}
private:
};
struct OpBoundaryResidualVector :
public OpEdgeEle {
public:
boost::shared_ptr<DataAtGaussPoints> &common_data)
if (nb_dofs) {
const double edge = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
double boundary_term =
for (int rr = 0; rr != nb_dofs; rr++) {
++t_base;
}
++t_w;
++t_coords;
++t_field;
}
ADD_VALUES);
}
}
private:
};
};
#endif
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEM::FaceElementForcesAndSourcesCore FaceEle
EntitiesFieldData::EntData EntData
#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.
EdgeElementForcesAndSourcesCoreSwitch< 0 > EdgeElementForcesAndSourcesCore
Edge finite element default.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
FTensor::Index< 'i', 2 > i
boost::function< double(const double, const double, const double)> ScalarFunc
constexpr auto VecSetValues
constexpr auto MatSetValues
const double r
rate factor
bool sYmm
If true assume that matrix is symmetric structure.
default operator for EDGE element
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.
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
default operator for TRI element
double getMeasure()
get measure of element
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
friend class UserDataOperator
VectorDouble locBoundaryRhs
OpBoundaryResidualVector(std::string field_name, ScalarFunc boundary_function, boost::shared_ptr< DataAtGaussPoints > &common_data)
boost::shared_ptr< DataAtGaussPoints > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MatrixDouble locBoundaryLhs
OpBoundaryTangentMatrix(std::string row_field_name, std::string col_field_name)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MatrixDouble transLocBoundaryLhs
boost::shared_ptr< DataAtGaussPoints > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
OpDomainResidualVector(std::string field_name, ScalarFunc source_term_function, boost::shared_ptr< DataAtGaussPoints > &common_data, boost::shared_ptr< std::vector< unsigned char >> boundary_marker=nullptr)
ScalarFunc sourceTermFunc
boost::shared_ptr< DataAtGaussPoints > commonData
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
OpDomainTangentMatrix(std::string row_field_name, std::string col_field_name, boost::shared_ptr< DataAtGaussPoints > &common_data, boost::shared_ptr< std::vector< unsigned char >> boundary_marker=nullptr)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Integrate the domain residual vector (RHS)
Integrate the domain tangent matrix (LHS)
#include <stdlib.h>
#include <BasicFiniteElements.hpp>
static char help[] =
"...\n\n";
public:
private:
static double sourceTermFunction(const double x, const double y,
const double z) {
return 200 * sin(x * 10.) * cos(y * 10.);
}
static double boundaryFunction(const double x, const double y,
const double z) {
return sin(x * 10.) * cos(y * 10.);
}
Simple *simpleInterface;
MPI_Comm mpiComm;
const int mpiRank;
SmartPetscObj<DM> dM;
SmartPetscObj<SNES> snesSolver;
std::string domainField;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
boost::shared_ptr<FaceEle> domainTangentMatrixPipeline;
boost::shared_ptr<FaceEle> domainResidualVectorPipeline;
boost::shared_ptr<EdgeEle> boundaryTangentMatrixPipeline;
boost::shared_ptr<EdgeEle> boundaryResidualVectorPipeline;
boost::shared_ptr<DataAtGaussPoints> previousUpdate;
boost::shared_ptr<VectorDouble> fieldValuePtr;
boost::shared_ptr<MatrixDouble> fieldGradPtr;
boost::shared_ptr<FaceEle> postProc;
Range boundaryEntitiesForFieldsplit;
};
: domainField("U"), mField(m_field), mpiComm(mField.get_comm()),
mpiRank(mField.get_comm_rank()) {
domainTangentMatrixPipeline = boost::shared_ptr<FaceEle>(
new FaceEle(mField));
domainResidualVectorPipeline =
boost::shared_ptr<FaceEle>(
new FaceEle(mField));
boundaryTangentMatrixPipeline =
boost::shared_ptr<EdgeEle>(
new EdgeEle(mField));
boundaryResidualVectorPipeline =
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);
}
}
}
}
auto domain_rule_lhs = [](
int,
int,
int p) ->
int {
return 2 * (
p + 1); };
auto domain_rule_rhs = [](
int,
int,
int p) ->
int {
return 2 * (
p + 1); };
auto boundary_rule_lhs = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
auto boundary_rule_rhs = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
}
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<unsigned char>>();
skin_edges, *marker_ptr);
return marker_ptr;
};
}
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
{
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
new OpSetInvJacH1ForFace(inv_jac_ptr));
}
{
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
new OpSetInvJacH1ForFace(inv_jac_ptr));
}
{
}
{
}
{
boost::shared_ptr<FaceEle> null_face;
null_face);
null_face);
boost::shared_ptr<EdgeEle> null_edge;
null_edge);
null_edge);
}
}
SmartPetscObj<Vec> global_rhs, global_solution;
if (1) {
KSP ksp_solver;
PC pc;
CHKERR KSPGetPC(ksp_solver, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
IS is_boundary;
CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
CHKERR ISDestroy(&is_boundary);
}
}
SCATTER_REVERSE);
}
CHKERR boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(
postProc)
->generateReferenceElementMesh();
CHKERR boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(
postProc)
CHKERR boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(
postProc)
->writeFile("out_result.h5m");
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
CHKERR poisson_problem.runProgram();
}
return 0;
}
int main(int argc, char *argv[])
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
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 _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.
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.
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.
keeps basic data about problem
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode outputResults()
MoFEMErrorCode boundaryCondition()
static double boundaryFunction(const double x, const double y, const double z)
MoFEMErrorCode solveSystem()
SmartPetscObj< SNES > snesSolver
static double sourceTermFunction(const double x, const double y, const double z)
MoFEMErrorCode readMesh()
MoFEM::Interface & mField
MoFEMErrorCode assembleSystem()
NonlinearPoisson(MoFEM::Interface &m_field)
boost::shared_ptr< EdgeEle > boundaryTangentMatrixPipeline
boost::shared_ptr< VectorDouble > fieldValuePtr
MoFEMErrorCode setIntegrationRules()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
boost::shared_ptr< MatrixDouble > fieldGradPtr
MoFEMErrorCode runProgram()
boost::shared_ptr< EdgeEle > boundaryResidualVectorPipeline
boost::shared_ptr< FaceEle > domainResidualVectorPipeline
MoFEMErrorCode setupProblem()
Range boundaryEntitiesForFieldsplit
boost::shared_ptr< DataAtGaussPoints > previousUpdate
boost::shared_ptr< FaceEle > postProc
boost::shared_ptr< FaceEle > domainTangentMatrixPipeline