12#include <boost/math/constants/constants.hpp>
14constexpr double third = boost::math::constants::third<double>();
19 OpHMHH(
const int tag,
const bool eval_rhs,
const bool eval_lhs,
20 boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
21 boost::shared_ptr<PhysicalEquations> &physics_ptr)
22 :
OpJacobian(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr) {}
37 const auto nb_integration_pts = getGaussPts().size2();
38 auto iu = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
39 auto t_grad_h1 = getFTensor2FromMat<3, 3>(
dataAtPts->wGradH1AtPts);
41 auto create_data_vec = [nb_integration_pts](
auto &
v) {
42 v.resize(nb_integration_pts,
false);
44 auto create_data_mat = [nb_integration_pts](
auto &
m) {
45 m.resize(9, nb_integration_pts,
false);
51 auto r_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
52 for (
unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
67 int r = ::function(
tAg,
physicsPtr->dependentVariablesPiola.size(),
73 "ADOL-C function evaluation with error r = %d", r);
106 const int number_of_active_variables =
physicsPtr->activeVariables.size();
107 const int number_of_dependent_variables =
109 std::vector<double *> jac_ptr(number_of_dependent_variables);
110 for (
unsigned int n = 0;
n != number_of_dependent_variables; ++
n) {
113 ->dependentVariablesPiolaDirevatives[
n *
114 number_of_active_variables]);
117 const auto nb_integration_pts = getGaussPts().size2();
119 auto create_data_mat = [nb_integration_pts](
auto &
m) {
120 m.resize(9, nb_integration_pts,
false);
123 dataAtPts->P_du.resize(81, nb_integration_pts,
false);
125 auto iu = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
126 auto t_grad_h1 = getFTensor2FromMat<3, 3>(
dataAtPts->wGradH1AtPts);
127 auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(
dataAtPts->P_du);
135 for (
unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
151 int r = ::jacobian(
tAg, number_of_dependent_variables,
152 number_of_active_variables,
153 &
physicsPtr->activeVariables[0], &jac_ptr[0]);
156 "ADOL-C function evaluation with error");
170 (t_P_dh_tmp(
i,
j,
o,
p) * t_h1_du(
o,
p,
m,
n)) * t_h1_du(
i,
j,
k,
l);
173 r_P_du(
i,
j,
m,
n) = t_P_dh_tmp(
i,
j,
m,
n);
199 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"stvenant_",
"",
"none");
201 CHKERR PetscOptionsScalar(
"-young_modulus",
"Young modulus",
"",
E, &
E,
203 CHKERR PetscOptionsScalar(
"-poisson_ratio",
"poisson ratio",
"", nu, &nu,
206 ierr = PetscOptionsEnd();
209 MOFEM_LOG(
"EP", Sev::inform) <<
"St Venant Kirchhoff model parameters: "
210 <<
"E = " <<
E <<
", nu = " << nu;
220 boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
221 boost::shared_ptr<PhysicalEquations> &physics_ptr) {
222 return (
new OpHMHH(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
267 ih(
i,
j) = (*t_h_ptr)(
i,
j);
274 enableMinMaxUsingAbs();
326 boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
327 boost::shared_ptr<PhysicalEquations> &physics_ptr) {
328 return (
new OpHMHH(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
333 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"mooneyrivlin_",
"",
"none");
336 CHKERR PetscOptionsScalar(
"-alpha",
"Alpha",
"",
alpha, &
alpha, PETSC_NULL);
338 CHKERR PetscOptionsScalar(
"-beta",
"Beta",
"",
beta, &
beta, PETSC_NULL);
348 ierr = PetscOptionsEnd();
404 ih(
i,
j) = (*t_h_ptr)(
i,
j);
411 enableMinMaxUsingAbs();
434 (levi_civita(
I,
J,
K) *
tF(
k,
K));
447 const int tape,
const double lambda,
const double mu,
448 const double sigma_y) {
456 const int tape,
const double alpha,
const double beta,
const double lambda,
457 const double sigma_y) {
460 boost::make_shared<HMHPMooneyRivlinWriggersEq63>(alpha, beta,
lambda);
474 boost::shared_ptr<HMHHencky> hencky_ptr)
477 "getOptions failed");
479 "Can not get data from block");
488 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
490 b.bulkModulusK, b.shearModulusG);
517 boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
518 boost::shared_ptr<PhysicalEquations> &physics_ptr) {
520 data_ptr, boost::dynamic_pointer_cast<HMHHencky>(physics_ptr)));
525 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"hencky_",
"",
"none");
527 CHKERR PetscOptionsScalar(
"-young_modulus",
"Young modulus",
"",
E, &
E,
529 CHKERR PetscOptionsScalar(
"-poisson_ratio",
"poisson ratio",
"",
nu, &
nu,
532 ierr = PetscOptionsEnd();
543 (boost::format(
"%s(.*)") %
"MAT_ELASTIC").str()
555 for (
auto m : meshset_vec_ptr) {
557 std::vector<double> block_data;
558 CHKERR m->getAttributes(block_data);
559 if (block_data.size() != 2) {
561 "Expected that block has two attribute");
563 auto get_block_ents = [&]() {
587 auto set_material_stiffness = [&]() {
597 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
605 set_material_stiffness();
Eshelbian plasticity interface.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#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_OPERATION_UNSUCCESSFUL
@ 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< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static constexpr auto size_symm
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr IntegrationType I
double young_modulus
Young modulus.
static enum RotSelector gradApperoximator
MoFEMErrorCode addMaterial_Hencky(double E, double nu)
MoFEMErrorCode addMaterial_HMHHStVenantKirchhoff(const int tape, const double lambda, const double mu, const double sigma_y)
MoFEMErrorCode addMaterial_HMHMooneyRivlin(const int tape, const double alpha, const double beta, const double lambda, const double sigma_y)
boost::shared_ptr< PhysicalEquations > physicalEquations
MoFEM::Interface & mField
MoFEMErrorCode evaluateRhs(EntData &data)
boost::shared_ptr< HMHHencky > henckyPtr
MoFEMErrorCode evaluateLhs(EntData &data)
OpHenckyJacobian(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< HMHHencky > hencky_ptr)
boost::shared_ptr< DataAtIntegrationPts > dataAtGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
virtual OpJacobian * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
MoFEMErrorCode extractBlockData(Sev sev)
std::vector< BlockData > blockData
MoFEMErrorCode extractBlockData(std::vector< const CubitMeshSets * > meshset_vec_ptr, Sev sev)
MoFEMErrorCode getMatDPtr(boost::shared_ptr< MatrixDouble > mat_D_ptr, double bulk_modulus_K, double shear_modulus_G)
HMHHencky(MoFEM::Interface &m_field, const double E, const double nu)
MoFEMErrorCode getOptions(boost::shared_ptr< DataAtIntegrationPts > data_ptr)
MoFEM::Interface & mField
MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h)
HMHPMooneyRivlinWriggersEq63(const double alpha, const double beta, const double lambda)
static constexpr int numberOfActiveVariables
virtual OpJacobian * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
static constexpr int numberOfDependentVariables
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
MoFEMErrorCode getOptions()
HMHStVenantKirchhoff(const double lambda, const double mu)
static constexpr int numberOfDependentVariables
virtual OpJacobian * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
MoFEMErrorCode getOptions()
static constexpr int numberOfActiveVariables
OpHMHH(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
MoFEMErrorCode evaluateLhs(EntData &data)
MoFEMErrorCode evaluateRhs(EntData &data)
boost::shared_ptr< PhysicalEquations > physicsPtr
material physical equations
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.