15#include <boost/math/constants/constants.hpp>
17constexpr double third = boost::math::constants::third<double>();
24 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
41 int nb_integration_pts = data.
getN().size1();
44 auto t_approx_P_adjoint_log_du =
56 auto get_ftensor2 = [](
auto &
v) {
58 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
61 int nb_base_functions = data.
getN().size2();
63 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
65 auto t_nf = get_ftensor2(
nF);
68 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
k,
l) * t_L(
k,
l,
L);
78 a * (t_approx_P_adjoint_log_du(
L) - t_Ldiff_u(
i,
j,
L) * t_P(
i,
j) -
79 t_L(
i,
j,
L) * t_viscous_P(
i,
j));
82 for (; bb != nb_dofs / 6; ++bb) {
83 t_nf(
L) -= t_row_base_fun * t_residual(
L);
87 for (; bb != nb_base_functions; ++bb)
91 ++t_approx_P_adjoint_log_du;
101 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha_u) {
107 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
108 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
109 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv) {
115 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
144 int nb_integration_pts = row_data.
getN().size1();
145 int row_nb_dofs = row_data.
getIndices().size();
146 int col_nb_dofs = col_data.
getIndices().size();
148 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
152 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
153 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
155 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
156 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
158 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
159 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
161 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
162 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
164 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
165 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
167 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
168 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
182 auto get_dP = [&]() {
188 auto t_approx_P_adjoint__dstretch =
200 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
203 t_deltaP(
i,
j) = t_approx_P_adjoint__dstretch(
i,
j) - t_P(
i,
j);
206 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
208 -t_Ldiff_u(
i,
j,
L) * (r_P_du(
i,
j,
k,
l) * t_Ldiff_u(
k,
l,
J));
211 (t_L(
i,
j,
L) * (t_diff(
i,
j,
k,
l) * t_L(
k,
l,
J)));
216 t_deltaP_sym(
i,
j) = (t_deltaP(
i,
j) || t_deltaP(
j,
i));
217 t_deltaP_sym(
i,
j) /= 2.0;
221 t_dP(
L,
J) += t_L(
i,
j,
L) * (t_diff2_uP2(
i,
j,
k,
l) * t_L(
k,
l,
J));
227 ++t_approx_P_adjoint__dstretch;
238 int row_nb_base_functions = row_data.
getN().size2();
241 auto t_dP = get_dP();
242 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
246 for (; rr != row_nb_dofs / 6; ++rr) {
248 auto t_m = get_ftensor2(
K, 6 * rr, 0);
249 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
250 const double b =
a * t_row_base_fun * t_col_base_fun;
251 t_m(
L,
J) -= b * t_dP(
L,
J);
258 for (; rr != row_nb_base_functions; ++rr) {
269 std::string row_field, std::string col_field,
270 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha) {
275 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
276 boost::shared_ptr<double> total_energy_ptr) {
281 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
282 boost::shared_ptr<PhysicalEquations> physics_ptr) {
287 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
288 boost::shared_ptr<PhysicalEquations> physics_ptr) {
309 boost::shared_ptr<double> scale_ptr,
310 boost::shared_ptr<PhysicalEquations> physics_ptr) {
315 OpHMHH(
const int tag,
const bool eval_rhs,
const bool eval_lhs,
316 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
317 boost::shared_ptr<PhysicalEquations> physics_ptr)
318 :
OpJacobian(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr) {}
333 const auto nb_integration_pts = getGaussPts().size2();
337 auto create_data_vec = [nb_integration_pts](
auto &
v) {
338 v.resize(nb_integration_pts,
false);
340 auto create_data_mat = [nb_integration_pts](
auto &
m) {
341 m.resize(9, nb_integration_pts,
false);
348 for (
unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
367 int r = ::function(
tAg,
physicsPtr->dependentVariablesPiola.size(),
373 "ADOL-C function evaluation with error r = %d", r);
409 const int number_of_active_variables =
physicsPtr->activeVariables.size();
410 const int number_of_dependent_variables =
412 std::vector<double *> jac_ptr(number_of_dependent_variables);
413 for (
unsigned int n = 0;
n != number_of_dependent_variables; ++
n) {
416 ->dependentVariablesPiolaDirevatives[
n *
417 number_of_active_variables]);
420 const auto nb_integration_pts = getGaussPts().size2();
422 auto create_data_mat = [nb_integration_pts](
auto &
m) {
423 m.resize(9, nb_integration_pts,
false);
426 dataAtPts->P_du.resize(81, nb_integration_pts,
false);
438 for (
unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
458 int r = ::jacobian(
tAg, number_of_dependent_variables,
459 number_of_active_variables,
460 &
physicsPtr->activeVariables[0], &jac_ptr[0]);
463 "ADOL-C function evaluation with error");
476 (t_P_dh_tmp(
i,
j, o, p) * t_h1_du(o, p,
m,
n)) * t_h1_du(
i,
j,
k,
l);
482 (t_P_dh_tmp(
i,
j, o, p) * t_h1_du(o, p,
m,
n)) * t_h1_du(
i,
j,
k,
l);
486 r_P_du(
i,
j,
m,
n) = t_P_dh_tmp(
i,
j,
m,
n);
498 const int tag,
const bool eval_rhs,
const bool eval_lhs,
499 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
500 boost::shared_ptr<PhysicalEquations> physics_ptr) {
501 return (
new OpHMHH(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
514 const int tape,
const double lambda,
const double mu,
515 const double sigma_y) {
523 const int tape,
const double alpha,
const double beta,
const double lambda,
524 const double sigma_y) {
527 boost::make_shared<HMHPMooneyRivlinWriggersEq63>(alpha, beta,
lambda);
Auxilary functions for Eshelbian plasticity.
Eshelbian plasticity interface.
Implementation StVenantKirchhoff.
Implementation of NeoHookean material.
Implement of MooneyRivlinWriggersEq63.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getDiffDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, C &&t_S, const int nb)
Get the Diff Diff Mat object.
static constexpr auto size_symm
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
static FTensor::Tensor4< FTensor::PackPtr< T *, 1 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3 > getFTensor4FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 4 (non symmetric) form data matrix.
constexpr auto field_name
FTensor::Index< 'm', 3 > m
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)
static enum StretchSelector stretchSelector
boost::shared_ptr< PhysicalEquations > physicalEquations
static enum RotSelector gradApproximator
static boost::function< double(const double)> f
static boost::function< double(const double)> dd_f
static boost::function< double(const double)> d_f
MoFEMErrorCode addMaterial_HMHNeohookean(const int tape, const double c10, const double K)
MoFEM::Interface & mField
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< double > scalePtr
OpGetScale(boost::shared_ptr< double > scale_ptr)
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
MoFEMErrorCode integratePiola(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
OpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
MoFEMErrorCode integrate(EntData &data)
virtual VolUserDataOperator * returnOpCalculateVarStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
virtual VolUserDataOperator * returnOpCalculateEnergy(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< double > total_energy_ptr)
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
virtual VolUserDataOperator * returnOpSpatialPhysicalExternalStrain(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< ExternalStrainVec > external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
virtual VolUserDataOperator * returnOpCalculateStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
virtual VolUserDataOperator * returnOpSetScale(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
virtual VolUserDataOperator * returnOpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
virtual UserDataOperator * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
bool sYmm
If true assume that matrix is symmetric structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
@ OPSPACE
operator do Work is execute on space data
double getVolume() const
element volume (linear geometry)