9#ifndef __KELVIN_VOIGT_DAMPER_HPP__
10#define __KELVIN_VOIGT_DAMPER_HPP__
13#error "MoFEM need to be compiled with ADOL-C"
56 MatrixBoundedArray<TYPE, 9>
F;
57 MatrixBoundedArray<TYPE, 9>
FDot;
58 MatrixBoundedArray<TYPE, 9>
61 MatrixBoundedArray<TYPE, 9>
63 MatrixBoundedArray<TYPE, 9>
65 MatrixBoundedArray<TYPE, 9>
invF;
91 for (
int ii = 0; ii < 3; ii++) {
119 for (
int ii = 0; ii < 3; ii++) {
139 invF.resize(3, 3,
false);
141 auto t_dashpotFirstPiolaKirchhoffStress = getFTensor2FromArray3by3(
143 auto t_dashpotCauchyStress = getFTensor2FromArray3by3(
153 J = determinantTensor3by3(t_F);
154 CHKERR invertTensor3by3(t_F,
J, t_invF);
155 t_dashpotFirstPiolaKirchhoffStress(
i,
j) =
156 J * (t_dashpotCauchyStress(
i,
k) * t_invF(
j,
k));
162 typedef boost::ptr_map<int, KelvinVoigtDamper::ConstitutiveEquation<adouble>>
263 bool calc_val,
bool calc_grad,
bool calc_dot =
false,
274 EntitiesFieldData::EntData &data) {
277 int nb_dofs = data.getFieldData().size();
281 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
282 int nb_gauss_pts = data.getN().size1();
290 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
297 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
320 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
321 for (
int rr1 = 0; rr1 < rank; rr1++) {
326 for (
int rr2 = 0; rr2 < 3; rr2++) {
328 t_diff_disp(rr1, rr2);
355 CommonData &common_data,
bool calculate_residual,
356 bool calculate_jacobian)
376 cE.
F.resize(3, 3,
false);
377 cE.
FDot.resize(3, 3,
false);
380 MatrixDouble &F_dot =
386 for (
int dd1 = 0; dd1 < 3; dd1++) {
387 for (
int dd2 = 0; dd2 < 3; dd2++) {
388 cE.
F(dd1, dd2) <<=
F(dd1, dd2);
392 for (
int dd1 = 0; dd1 < 3; dd1++) {
393 for (
int dd2 = 0; dd2 < 3; dd2++) {
394 cE.
FDot(dd1, dd2) <<= F_dot(dd1, dd2);
408 for (
int d1 = 0; d1 < 3; d1++) {
409 for (
int d2 = 0; d2 < 3; d2++) {
429 "ADOL-C function evaluation with error r = %d", r);
445 "ADOL-C function evaluation with error");
447 }
catch (
const std::exception &ex) {
448 std::ostringstream ss;
449 ss <<
"throw in method: " << ex.what() << std::endl;
450 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
467 MatrixDouble &F_dot =
469 int nb_active_variables = 0;
472 for (
int dd1 = 0; dd1 < 3; dd1++) {
473 for (
int dd2 = 0; dd2 < 3; dd2++) {
478 for (
int dd1 = 0; dd1 < 3; dd1++) {
479 for (
int dd2 = 0; dd2 < 3; dd2++) {
486 "Number of active variables does not much");
518 EntitiesFieldData::EntData &row_data) {
521 if (row_type != MBVERTEX)
549 EntitiesFieldData::EntData &row_data) {
551 int nb_dofs = row_data.getIndices().size();
552 int *indices_ptr = &row_data.getIndices()[0];
569 EntitiesFieldData::EntData &row_data) {
576 int nb_dofs = row_data.getIndices().size();
580 nF.resize(nb_dofs,
false);
582 int nb_gauss_pts = row_data.getN().size1();
583 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
584 const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_dofs / 3);
585 const MatrixDouble &stress =
588 for (
int dd = 0; dd < nb_dofs / 3; dd++) {
589 for (
int rr = 0; rr < 3; rr++) {
590 for (
int nn = 0; nn < 3; nn++) {
591 nF[3 * dd + rr] += val * diffN(dd, nn) * stress(rr, nn);
610 EntitiesFieldData::EntData &row_data,
611 EntitiesFieldData::EntData &col_data) {
613 int nb_row = row_data.getIndices().size();
614 int nb_col = col_data.getIndices().size();
615 int *row_indices_ptr = &row_data.getIndices()[0];
616 int *col_indices_ptr = &col_data.getIndices()[0];
618 col_indices_ptr, &
K(0, 0), ADD_VALUES);
621 if (row_side != col_side || row_type != col_type) {
622 transK.resize(nb_col, nb_row,
false);
625 nb_row, row_indices_ptr, &
transK(0, 0),
639 common_data.spatialPositionName),
645 int nb_col = col_data.getIndices().size();
648 const MatrixAdaptor diffN = col_data.getDiffN(gg, nb_col / 3);
650 for (
int dd = 0; dd < nb_col / 3; dd++) {
651 for (
int jj = 0; jj < 3; jj++) {
652 double a = diffN(dd, jj);
653 for (
int rr = 0; rr < 3; rr++) {
654 for (
int ii = 0; ii < 9;
656 dStress_dx(ii, 3 * dd + rr) += jac_stress(ii, 3 * rr + jj) *
a;
665 EntitiesFieldData::EntData &row_data,
666 EntitiesFieldData::EntData &col_data) {
673 int nb_row = row_data.getIndices().size();
674 int nb_col = col_data.getIndices().size();
679 K.resize(nb_row, nb_col,
false);
681 int nb_gauss_pts = row_data.getN().size1();
682 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
687 const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_row / 3);
689 for (
int dd1 = 0; dd1 < nb_row / 3; dd1++) {
690 for (
int rr1 = 0; rr1 < 3; rr1++) {
691 for (
int dd2 = 0; dd2 < nb_col / 3; dd2++) {
692 for (
int rr2 = 0; rr2 < 3; rr2++) {
693 K(3 * dd1 + rr1, 3 * dd2 + rr2) +=
694 (diffN(dd1, 0) *
dStress_dx(3 * rr1 + 0, 3 * dd2 + rr2) +
695 diffN(dd1, 1) *
dStress_dx(3 * rr1 + 1, 3 * dd2 + rr2) +
696 diffN(dd1, 2) *
dStress_dx(3 * rr1 + 2, 3 * dd2 + rr2));
705 CHKERR aSemble(row_side, col_side, row_type, col_type, row_data,
718 common_data.spatialPositionName),
724 int nb_col = col_data.getIndices().size();
727 const MatrixAdaptor diffN = col_data.getDiffN(gg, nb_col / 3);
729 for (
int dd = 0; dd < nb_col / 3; dd++) {
730 for (
int jj = 0; jj < 3; jj++) {
731 double a = diffN(dd, jj);
732 for (
int rr = 0; rr < 3; rr++) {
733 for (
int ii = 0; ii < 9;
745 EntitiesFieldData::EntData &row_data,
746 EntitiesFieldData::EntData &col_data) {
753 int nb_row = row_data.getIndices().size();
754 int nb_col = col_data.getIndices().size();
759 K.resize(nb_row, nb_col,
false);
761 int nb_gauss_pts = row_data.getN().size1();
762 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
767 const MatrixAdaptor &diffN = row_data.getDiffN(gg, nb_row / 3);
769 for (
int dd1 = 0; dd1 < nb_row / 3; dd1++) {
770 for (
int rr1 = 0; rr1 < 3; rr1++) {
771 for (
int dd2 = 0; dd2 < nb_col / 3; dd2++) {
772 for (
int rr2 = 0; rr2 < 3; rr2++) {
773 K(3 * dd1 + rr1, 3 * dd2 + rr2) +=
774 (diffN(dd1, 0) *
dStress_dot(3 * rr1 + 0, 3 * dd2 + rr2) +
775 diffN(dd1, 1) *
dStress_dot(3 * rr1 + 1, 3 * dd2 + rr2) +
776 diffN(dd1, 2) *
dStress_dot(3 * rr1 + 2, 3 * dd2 + rr2));
785 CHKERR aSemble(row_side, col_side, row_type, col_type, row_data,
796 if (it->getName().compare(0, 6,
"DAMPER") == 0) {
797 std::vector<double> data;
798 CHKERR it->getAttributes(data);
799 if (data.size() < 2) {
800 SETERRQ(PETSC_COMM_SELF, 1,
"Data inconsistency");
818 fe_ptr->getOpPtrVector().push_back(
824 fe_ptr->getOpPtrVector().push_back(
825 new OpCalculateVectorFieldGradientDot<3, 3>(
833 std::vector<int> tags;
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#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 ...
@ MOFEM_OPERATION_UNSUCCESSFUL
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#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
constexpr auto field_name
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
AssembleMatrix(string row_name, string col_name)
MoFEMErrorCode aSemble(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
AssembleVector(string field_name)
Dumper material parameters.
double vBeta
Poisson ration spring alpha.
double gBeta
Sheer modulus spring alpha.
Common data for nonlinear_elastic_elem model.
std::map< int, int > nbActiveVariables
boost::shared_ptr< MatrixDouble > dataAtGaussTmpPtr
std::vector< MatrixDouble > dashpotFirstPiolaKirchhoffStress
std::vector< double * > jacRowPtr
string spatialPositionNameDot
std::vector< MatrixDouble > jacStress
string spatialPositionName
string meshNodePositionName
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::map< int, int > nbActiveResults
boost::shared_ptr< MatrixDouble > gradDataAtGaussTmpPtr
Constitutive model functions.
FTensor::Index< 'j', 3 > j
virtual ~ConstitutiveEquation()=default
MatrixBoundedArray< TYPE, 9 > FDot
Rate of gradient of deformation.
FTensor::Index< 'k', 3 > k
MatrixBoundedArray< TYPE, 9 > invF
Inverse of gradient of deformation.
FTensor::Index< 'i', 3 > i
MatrixBoundedArray< TYPE, 9 > F
Gradient of deformation.
virtual MoFEMErrorCode calculateFirstPiolaKirchhoffStress()
Calculate First Piola-Kirchhoff Stress Dashpot stress.
MatrixBoundedArray< TYPE, 9 > gradientUDot
Rate of gradient of displacements.
virtual MoFEMErrorCode calculateDashpotCauchyStress()
Calculate Cauchy dashpot stress.
MatrixBoundedArray< TYPE, 9 > engineringStrainDot
virtual MoFEMErrorCode calculateEngineeringStrainDot()
Calculate strain rate.
MatrixBoundedArray< TYPE, 9 > dashpotFirstPiolaKirchhoffStress
Stress generated by spring beta.
ConstitutiveEquation(BlockMaterialData &data, bool is_displacement=true)
MatrixBoundedArray< TYPE, 9 > dashpotCauchyStress
Stress generated by spring beta.
TYPE traceEngineeringStrainDot
TYPE J
Jacobian of gradient of deformation.
definition of volume element
int addToRule
Takes into account HO geometry.
MoFEMErrorCode postProcess()
function is run at the end of loop
DamperFE(MoFEM::Interface &m_field, CommonData &common_data)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
OpGetDataAtGaussPts(const std::string field_name, CommonData &common_data, bool calc_val, bool calc_grad, bool calc_dot=false, EntityType zero_at_type=MBVERTEX)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator field value.
std::map< int, int > & nbActiveVariables
MoFEMErrorCode calculateJacobian(TagEvaluate te)
bool calculateJacobianBool
bool calculateResidualBool
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
std::map< int, int > & nbActiveResults
KelvinVoigtDamper::ConstitutiveEquation< adouble > & cE
MoFEMErrorCode calculateAtIntPtsDamperStress()
VectorDouble activeVariables
MoFEMErrorCode calculateFunction(TagEvaluate te, double *ptr)
OpJacobian(const std::string field_name, std::vector< int > tags, KelvinVoigtDamper::ConstitutiveEquation< adouble > &ce, CommonData &common_data, bool calculate_residual, bool calculate_jacobian)
MoFEMErrorCode recordDamperStress()
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode get_dStress_dot(EntitiesFieldData::EntData &col_data, int gg)
OpLhsdxdot(CommonData &common_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpLhsdxdx(CommonData &common_data)
MoFEMErrorCode get_dStress_dx(EntitiesFieldData::EntData &col_data, int gg)
Assemble internal force vector.
OpRhsStress(CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Implementation of Kelvin Voigt Damper.
MoFEMErrorCode setBlockDataMap()
ConstitutiveEquationMap constitutiveEquationMap
std::map< int, BlockMaterialData > blockMaterialDataMap
KelvinVoigtDamper(MoFEM::Interface &m_field)
MoFEMErrorCode setOperators(const int tag)
boost::ptr_map< int, KelvinVoigtDamper::ConstitutiveEquation< adouble > > ConstitutiveEquationMap
MoFEM::Interface & mField
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
@ OPCOL
operator doWork function is executed on FE columns
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
double getVolume() const
element volume (linear geometry)
Volume finite element base.
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)