70 "extract block data failed");
72#ifdef NEOHOOKEAN_SCALING
87 "Stretch selector is not equal to LOG");
91 "Exponent base is not equal to exp(1)");
99 if (b.blockEnts.find(ent) != b.blockEnts.end()) {
100 return std::make_pair(b.c10, b.K);
105 "Block not found for entity handle. If you mat set "
106 "block, set it to all elements");
112 boost::shared_ptr<PhysicalEquations> physics_ptr)
116 MoFEMErrorCode
doWork(
int side, EntityType type,
117 EntitiesFieldData::EntData &data) {
119 auto neohookean_ptr =
120 boost::dynamic_pointer_cast<HMHNeohookean>(
physicsPtr);
121 *(
scalePtr) = neohookean_ptr->eqScaling;
132 boost::shared_ptr<PhysicalEquations> physics_ptr) {
134 return new OpGetScale(scale_ptr, physics_ptr);
145 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
146 boost::shared_ptr<PhysicalEquations> physics_ptr) {
149 data_ptr, physics_ptr);
151 return (
new OpJacobian(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
157 PetscOptionsBegin(PETSC_COMM_WORLD,
"neo_hookean_",
"",
"none");
167 CHKERR PetscOptionsScalar(
"-viscosity_alpha_grad_u",
"viscosity",
"",
191 (boost::format(
"%s(.*)") %
"MAT_NEOHOOKEAN").str()
203 for (
auto m : meshset_vec_ptr) {
205 std::vector<double> block_data;
206 CHKERR m->getAttributes(block_data);
207 if (block_data.size() < 2) {
209 "Expected that block has atleast two attributes");
211 auto get_block_ents = [&]() {
217 double c10 = block_data[0];
218 double K = block_data[1];
220 blockData.push_back({c10, K, get_block_ents()});
222 MOFEM_LOG(
"EP", sev) <<
"MatBlock Neo-Hookean c10 = "
224 <<
" K = " <<
blockData.back().K <<
" nb ents. = "
245 ih(
i,
j) = (*t_h_ptr)(
i,
j);
252 enableMinMaxUsingAbs();
283 "gradApproximator not handled");
292 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
293 const double alpha_u);
303 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
304 const double alpha_u) {
317 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
318 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
319 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
330 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
331 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
332 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv) {
334 external_strain_vec_ptr, smv);
339 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
348 std::string row_field, std::string col_field,
349 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha) {
352 row_field, col_field, data_ptr, alpha);
386 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha_u)
391 double circumcenter[3],
double *xi,
double *
eta,
398 auto neohookean_ptr =
399 boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
400 if (!neohookean_ptr) {
402 "Pointer to HMHNeohookean is null");
404 auto [def_c10, def_K] =
405 neohookean_ptr->getMaterialParameters(getFEEntityHandle());
407 double c10 = def_c10 / neohookean_ptr->eqScaling;
408 double alpha_u = alphaU / neohookean_ptr->eqScaling;
409 double K = def_K / neohookean_ptr->eqScaling;
411 auto time_step = getTStimeStep();
412 double alpha_grad_u =
413 neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;
422 int nb_integration_pts = data.
getN().size1();
423 auto v = getVolume();
424 auto t_w = getFTensor0IntegrationWeight();
425 auto t_approx_P_adjoint_log_du =
426 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
427 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
428 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
430 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
432 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
434 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
436 getFTensor2FromMat<6, 3>(dataAtPts->gradLogStretchDotTensorAtPts);
438 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
440 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
441 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
442 auto &nbUniq = dataAtPts->nbUniq;
455 auto get_ftensor2 = [](
auto &
v) {
457 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
460 int nb_base_functions = data.
getN().size2();
467 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
472 auto t_neohookean_hencky =
478 const double tr = t_log_u(
i,
i);
480 t_P(L) = t_L(
i,
j, L) * (t_neohookean_hencky(
i,
j) +
483 t_viscous_P(L) = alpha_u * (t_L(
i,
j, L) * t_dot_log_u(
i,
j));
486 t_residual(L) = t_approx_P_adjoint_log_du(L) - t_P(L) - t_viscous_P(L);
490 t_grad_residual(L,
i) = alpha_grad_u * t_grad_log_u(L,
i);
491 t_grad_residual(L,
i) *=
a;
493 ++t_approx_P_adjoint_log_du;
498 auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
500 for (; bb != nb_dofs /
size_symm; ++bb) {
501 t_nf(L) -= t_row_base_fun * t_residual(L);
502 t_nf(L) += t_grad_base_fun(
i) * t_grad_residual(L,
i);
507 for (; bb != nb_base_functions; ++bb) {
519 "Not implemented for Neo-Hookean (used ADOL-C)");
533 "gradApproximator not handled");
541 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
542 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
543 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
545 externalStrainVecPtr(external_strain_vec_ptr), scalingMethodsMap{smv} {}
551 auto neohookean_ptr =
552 boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
553 if (!neohookean_ptr) {
555 "Pointer to HMHNeohookean is null");
566 for (
auto &ext_strain_block : (*externalStrainVecPtr)) {
568 if (ext_strain_block.ents.find(fe_ent) != ext_strain_block.ents.end()) {
571 if (scalingMethodsMap.find(ext_strain_block.blockName) !=
572 scalingMethodsMap.end()) {
574 scalingMethodsMap.at(ext_strain_block.blockName)->getScale(time);
577 <<
"No scaling method found for " << ext_strain_block.blockName;
580 auto [def_c10, def_K] =
581 neohookean_ptr->getMaterialParameters(getFEEntityHandle());
582 double c10 = def_c10 / neohookean_ptr->eqScaling;
584 double external_strain_val =
scale * ext_strain_block.val;
585 double K = ext_strain_block.bulkModulusK;
592 int nb_integration_pts = data.
getN().size1();
593 auto v = getVolume();
594 auto t_w = getFTensor0IntegrationWeight();
603 int nb_base_functions = data.
getN().size2();
606 const double tr = 3.0 * external_strain_val;
607 const double sigma_J = K * tr;
609 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
615 t_residual(L) += (t_L(
i,
j, L) *
t_kd(
i,
j)) * sigma_J;
618 auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
620 for (; bb != nb_dofs /
size_symm; ++bb) {
621 t_nf(L) -= t_row_base_fun * t_residual(L);
625 for (; bb != nb_base_functions; ++bb) {
635 std::string row_field, std::string col_field,
636 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha)
648 auto neohookean_ptr =
649 boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
650 if (!neohookean_ptr) {
652 "Pointer to HMHNeohookean is null");
654 auto [def_c10, def_K] =
655 neohookean_ptr->getMaterialParameters(getFEEntityHandle());
657 double c10 = def_c10 / neohookean_ptr->eqScaling;
658 double alpha_u = alphaU / neohookean_ptr->eqScaling;
659 double lambda = def_K / neohookean_ptr->eqScaling;
661 auto time_step = getTStimeStep();
662 double alpha_grad_u =
663 neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;
674 int nb_integration_pts = row_data.
getN().size1();
675 int row_nb_dofs = row_data.
getIndices().size();
676 int col_nb_dofs = col_data.
getIndices().size();
678 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
682 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
683 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
685 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
686 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
688 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
689 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
691 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
692 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
694 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
695 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
697 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
698 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
710 auto v = getVolume();
711 auto ts_a = getTSa();
712 auto t_w = getFTensor0IntegrationWeight();
714 int row_nb_base_functions = row_data.
getN().size2();
718 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
720 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
722 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
724 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
725 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
726 auto t_approx_P_adjoint__dstretch =
727 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
728 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
729 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
730 auto &nbUniq = dataAtPts->nbUniq;
737 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
742 auto d_neohookean = [c10,
lambda](
auto v) {
747 t_eigen_vals, t_eigen_vecs, neohookean, d_neohookean, t_nb_uniq);
749 const auto tr = t_log_u(
i,
i);
751 t_dP(L,
J) = -t_L(
i,
j, L) * ((t_diff_neohookean(
i,
j,
k,
l) +
753 t_kd_sym(
i,
j) * t_kd_sym(
k,
l)) *
755 t_dP(L,
J) -= (alpha_u * ts_a) *
756 (t_L(
i,
j, L) * (t_diff(
i,
j,
k,
l) * t_L(
k,
l,
J)));
760 t_deltaP(
i,
j) = (t_approx_P_adjoint__dstretch(
i,
j) ||
761 t_approx_P_adjoint__dstretch(
j,
i)) /
766 t_dP(L,
J) += t_L(
i,
j, L) * (t_diff2_uP(
i,
j,
k,
l) * t_L(
k,
l,
J));
768 ++t_approx_P_adjoint__dstretch;
775 for (; rr != row_nb_dofs /
size_symm; ++rr) {
779 auto t_m = get_ftensor2(K, 6 * rr, 0);
780 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
781 double b =
a * t_row_base_fun * t_col_base_fun;
782 double c = (
a * alpha_grad_u * ts_a) *
783 (t_row_grad_fun(
i) * t_col_grad_fun(
i));
784 t_m(L,
J) -= b * t_dP(L,
J);
785 t_m(L,
J) +=
c * t_kd_sym(L,
J);
795 for (; rr != row_nb_base_functions; ++rr) {
807 "Not implemented for Neo-Hookean (used ADOL-C)");
821 "gradApproximator not handled");
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define FTENSOR_INDEX(DIM, I)
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
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 getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
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.
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
static constexpr auto size_symm
constexpr IntegrationType I
double zeta
Viscous hardening.
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static double dynamicTime
static enum StretchSelector stretchSelector
static double exponentBase
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
static PetscBool dynamicRelaxation
OpGetScale(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
boost::shared_ptr< double > scalePtr
boost::shared_ptr< PhysicalEquations > physicsPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode evaluateRhs(EntData &data)
MoFEMErrorCode evaluateLhs(EntData &data)
MoFEMErrorCode integrate(EntData &data)
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
OpSpatialPhysicalExternalStrain(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)
OpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
MoFEMErrorCode integrate(EntData &data)
static double fun_diff_neohookean_bulk(double K, double tr)
Definition of derivative of axiator of Neo-hookean function.
static double fun_d_neohookean(double c10, double v)
Definition of derivative of Neo-hookean function.
static double fun_neohookean_bulk(double K, double tr)
Definition of axiator of Neo-hookean function.
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
std::vector< BlockData > blockData
MoFEMErrorCode extractBlockData(std::vector< const CubitMeshSets * > meshset_vec_ptr, Sev sev)
MoFEMErrorCode extractBlockData(Sev sev)
static constexpr int numberOfDependentVariables
VolUserDataOperator * returnOpSetScale(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
static double fun_neohookean(double c10, double v)
Definition of Neo-hookean function.
MoFEMErrorCode getOptions()
MoFEM::Interface & mField
auto getMaterialParameters(EntityHandle ent)
static constexpr int numberOfActiveVariables
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)
HMHNeohookean(MoFEM::Interface &m_field, const double c10, const double K)
VolUserDataOperator * returnOpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
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)
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
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)
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
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.
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.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.