v0.16.3
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | List of all members
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation Struct Reference

Constitutive (physical) equation for interface. More...

#include "tutorials/cor-12_cohesive_interface/src/CohesiveInterfaceElement.hpp"

Collaboration diagram for CohesiveElement::CohesiveInterfaceElement::PhysicalEquation:
[legend]

Public Member Functions

 PhysicalEquation (MoFEM::Interface &m_field)
 
virtual ~PhysicalEquation ()
 
MoFEMErrorCode iNitailise (const FEMethod *fe_method)
 Initialize history variable data.
 
double calcG (int gg, MatrixDouble gap_loc)
 Calculate gap opening.
 
MoFEMErrorCode getKappa (int nb_gauss_pts, const FEMethod *fe_method)
 Get pointer from the mesh to histoy variables \(\kappa\).
 
MoFEMErrorCode calcDglob (const double omega, MatrixDouble &R)
 Calculate stiffness material matrix.
 
MoFEMErrorCode calcOmega (const double kappa, double &omega)
 Calculate damage.
 
MoFEMErrorCode calcTangetDglob (const double omega, double g, const VectorDouble &gap_loc, MatrixDouble &R)
 Calculate tangent material stiffness.
 
virtual MoFEMErrorCode calculateTraction (VectorDouble &traction, int gg, CommonData &common_data, const FEMethod *fe_method)
 Calculate tractions.
 
virtual MoFEMErrorCode calculateTangentStiffeness (MatrixDouble &tangent_matrix, int gg, CommonData &common_data, const FEMethod *fe_method)
 Calculate tangent stiffness.
 
virtual MoFEMErrorCode updateHistory (CommonData &common_data, const FEMethod *fe_method)
 Update history variables when converged.
 

Public Attributes

MoFEM::Interface & mField
 
bool isInitialised
 
double h
 
double youngModulus
 
double beta
 
double ft
 
double Gf
 
Range pRisms
 
Tag thKappa
 
Tag thDamagedPrism
 
double E0
 
double g0
 
double kappa1
 
double * kappaPtr
 
int kappaSize
 
MatrixDouble Dglob
 
MatrixDouble Dloc
 

Detailed Description

Constitutive (physical) equation for interface.

This is linear degradation model. Material parameters are: strength \(f_t\), interface fracture energy \(G_f\), elastic material stiffness \(E\). Parameter \(\beta\) controls how interface opening is calculated.

Model parameter is interface penalty thickness \(h\).

Examples
mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 50 of file CohesiveInterfaceElement.hpp.

Constructor & Destructor Documentation

◆ PhysicalEquation()

CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::PhysicalEquation ( MoFEM::Interface &  m_field)
inline

Definition at line 54 of file CohesiveInterfaceElement.hpp.

55 : mField(m_field), isInitialised(false){};

◆ ~PhysicalEquation()

virtual CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::~PhysicalEquation ( )
inlinevirtual

Definition at line 57 of file CohesiveInterfaceElement.hpp.

57{}

Member Function Documentation

◆ calcDglob()

MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcDglob ( const double  omega,
MatrixDouble &  R 
)
inline

Calculate stiffness material matrix.

\[ \mathbf{D}_\textrm{loc} = (1-\Omega) \mathbf{I} E_0 \]

where \(E_0\) is initial interface penalty stiffness

\[ \mathbf{D}_\textrm{glob} = \mathbf{R}^\textrm{T} \mathbf{D}_\textrm{loc}\mathbf{R} \]

Definition at line 139 of file CohesiveInterfaceElement.hpp.

139 {
140 MoFEMFunctionBeginHot;
141 Dglob.resize(3, 3);
142 Dloc.resize(3, 3);
143 Dloc.clear();
144 double E = (1 - omega) * E0;
145 Dloc(0, 0) = E;
146 Dloc(1, 1) = E;
147 Dloc(2, 2) = E;
148 Dglob = prod(Dloc, R);
149 Dglob = prod(trans(R), Dglob);
150 MoFEMFunctionReturnHot(0);
151 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
omega
constexpr double omega
Save field DOFS on vertices/tags.
Definition dynamic_first_order_con_law.cpp:92
R
@ R
Definition free_surface.cpp:630

◆ calcG()

double CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcG ( int  gg,
MatrixDouble  gap_loc 
)
inline

Calculate gap opening.

\[ g = \sqrt{ g_n^2 + \beta(g_{s1}^2 + g_{s2}^2)} \]

Definition at line 94 of file CohesiveInterfaceElement.hpp.

94 {
95 return sqrt(pow(gap_loc(gg, 0), 2) +
96 beta * (pow(gap_loc(gg, 1), 2) + pow(gap_loc(gg, 2), 2)));
97 }

◆ calcOmega()

MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcOmega ( const double  kappa,
double &  omega 
)
inline

Calculate damage.

\[ \Omega = \frac{1}{2} \frac{(2 G_f E_0+f_t^2)\kappa}{(ft+E_0 \kappa)G_f} \]

Definition at line 160 of file CohesiveInterfaceElement.hpp.

160 {
161 MoFEMFunctionBeginHot;
162 omega = 0;
163 if (kappa >= kappa1) {
164 omega = 1;
165 MoFEMFunctionReturnHot(0);
166 } else if (kappa > 0) {
167 double a = (2.0 * Gf * E0 + ft * ft) * kappa;
168 double b = (ft + E0 * kappa) * Gf;
169 omega = 0.5 * a / b;
170 }
171 MoFEMFunctionReturnHot(0);
172 }
a
constexpr double a
Definition approx_sphere.cpp:30
kappa
double kappa
Definition free_surface.cpp:205

◆ calcTangetDglob()

MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcTangetDglob ( const double  omega,
double  g,
const VectorDouble &  gap_loc,
MatrixDouble &  R 
)
inline

Calculate tangent material stiffness.

Definition at line 176 of file CohesiveInterfaceElement.hpp.

178 {
179 MoFEMFunctionBeginHot;
180 Dglob.resize(3, 3);
181 Dloc.resize(3, 3);
182 double domega =
183 0.5 * (2 * Gf * E0 + ft * ft) / ((ft + (g - ft / E0) * E0) * Gf) -
184 0.5 * ((g - ft / E0) * (2 * Gf * E0 + ft * ft) * E0) /
185 (pow(ft + (g - ft / E0) * E0, 2) * Gf);
186 Dloc.resize(3, 3);
187 //r0
188 Dloc(0, 0) = (1 - omega) * E0 - domega * E0 * gap_loc[0] * gap_loc[0] / g;
189 Dloc(0, 1) = -domega * E0 * gap_loc[0] * beta * gap_loc[1] / g;
190 Dloc(0, 2) = -domega * E0 * gap_loc[0] * beta * gap_loc[2] / g;
191 //r1
192 Dloc(1, 0) = -domega * E0 * gap_loc[1] * gap_loc[0] / g;
193 Dloc(1, 1) =
194 (1 - omega) * E0 - domega * E0 * gap_loc[1] * beta * gap_loc[1] / g;
195 Dloc(1, 2) = -domega * E0 * gap_loc[1] * beta * gap_loc[2] / g;
196 //r2
197 Dloc(2, 0) = -domega * E0 * gap_loc[2] * gap_loc[0] / g;
198 Dloc(2, 1) = -domega * E0 * gap_loc[2] * beta * gap_loc[1] / g;
199 Dloc(2, 2) =
200 (1 - omega) * E0 - domega * E0 * gap_loc[2] * beta * gap_loc[2] / g;
201 Dglob = prod(Dloc, R);
202 Dglob = prod(trans(R), Dglob);
203 MoFEMFunctionReturnHot(0);
204 }
g
constexpr double g
Definition shallow_wave.cpp:63

◆ calculateTangentStiffeness()

virtual MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calculateTangentStiffeness ( MatrixDouble &  tangent_matrix,
int  gg,
CommonData &  common_data,
const FEMethod *  fe_method 
)
inlinevirtual

Calculate tangent stiffness.

Definition at line 239 of file CohesiveInterfaceElement.hpp.

241 {
242 MoFEMFunctionBeginHot;
243
244 try {
245 if (!isInitialised) {
246 CHKERR iNitailise(fe_method);
247 isInitialised = true;
248 }
249 if (gg == 0) {
250 CHKERR getKappa(common_data.gapGlob.size1(), fe_method);
251 }
252 double g = calcG(gg, common_data.gapLoc);
253 double kappa = fmax(g - g0, kappaPtr[gg]);
254 double omega = 0;
255 CHKERR calcOmega(kappa, omega);
256 int iter;
257 CHKERR SNESGetIterationNumber(fe_method->snes, &iter);
258 if ((kappa <= kappaPtr[gg]) || (kappa >= kappa1) || (iter <= 1)) {
259 CHKERR calcDglob(omega, common_data.R[gg]);
260 } else {
261 ublas::matrix_row<MatrixDouble> g_loc(common_data.gapLoc, gg);
262 CHKERR calcTangetDglob(omega, g, g_loc, common_data.R[gg]);
263 }
264 tangent_matrix.resize(3, 3);
265 noalias(tangent_matrix) = Dglob;
266 } catch (const std::exception &ex) {
267 std::ostringstream ss;
268 ss << "throw in method: " << ex.what() << std::endl;
269 SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, "%s",
270 ss.str().c_str());
271 }
272 MoFEMFunctionReturnHot(0);
273 }
MOFEM_STD_EXCEPTION_THROW
@ MOFEM_STD_EXCEPTION_THROW
Definition definitions.h:39
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcG
double calcG(int gg, MatrixDouble gap_loc)
Calculate gap opening.
Definition CohesiveInterfaceElement.hpp:94
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::iNitailise
MoFEMErrorCode iNitailise(const FEMethod *fe_method)
Initialize history variable data.
Definition CohesiveInterfaceElement.hpp:70
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcOmega
MoFEMErrorCode calcOmega(const double kappa, double &omega)
Calculate damage.
Definition CohesiveInterfaceElement.hpp:160
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcTangetDglob
MoFEMErrorCode calcTangetDglob(const double omega, double g, const VectorDouble &gap_loc, MatrixDouble &R)
Calculate tangent material stiffness.
Definition CohesiveInterfaceElement.hpp:176
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcDglob
MoFEMErrorCode calcDglob(const double omega, MatrixDouble &R)
Calculate stiffness material matrix.
Definition CohesiveInterfaceElement.hpp:139
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::getKappa
MoFEMErrorCode getKappa(int nb_gauss_pts, const FEMethod *fe_method)
Get pointer from the mesh to histoy variables .
Definition CohesiveInterfaceElement.hpp:104

◆ calculateTraction()

virtual MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calculateTraction ( VectorDouble &  traction,
int  gg,
CommonData &  common_data,
const FEMethod *  fe_method 
)
inlinevirtual

Calculate tractions.

\[ \mathbf{t} = \mathbf{D}_\textrm{glob}\mathbf{g} \]

Definition at line 213 of file CohesiveInterfaceElement.hpp.

215 {
216 MoFEMFunctionBeginHot;
217
218 if (!isInitialised) {
219 CHKERR iNitailise(fe_method);
220 isInitialised = true;
221 }
222 if (gg == 0) {
223 CHKERR getKappa(common_data.gapGlob.size1(), fe_method);
224 }
225 double g = calcG(gg, common_data.gapLoc);
226 double kappa = fmax(g - g0, kappaPtr[gg]);
227 double omega = 0;
228 CHKERR calcOmega(kappa, omega);
229 CHKERR calcDglob(omega, common_data.R[gg]);
230 traction.resize(3);
231 ublas::matrix_row<MatrixDouble> gap_glob(common_data.gapGlob, gg);
232 noalias(traction) = prod(Dglob, gap_glob);
233 MoFEMFunctionReturnHot(0);
234 }

◆ getKappa()

MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::getKappa ( int  nb_gauss_pts,
const FEMethod *  fe_method 
)
inline

Get pointer from the mesh to histoy variables \(\kappa\).

Definition at line 104 of file CohesiveInterfaceElement.hpp.

104 {
105 MoFEMFunctionBeginHot;
106 EntityHandle ent = fe_method->numeredEntFiniteElementPtr->getEnt();
107
108 rval = mField.get_moab().tag_get_by_ptr(
109 thKappa, &ent, 1, (const void **)&kappaPtr, &kappaSize);
110 if (rval != MB_SUCCESS || kappaSize != nb_gauss_pts) {
111 VectorDouble kappa;
112 kappa.resize(nb_gauss_pts);
113 kappa.clear();
114 int tag_size[1];
115 tag_size[0] = nb_gauss_pts;
116 void const *tag_data[] = {&kappa[0]};
117 CHKERR mField.get_moab().tag_set_by_ptr(thKappa, &ent, 1, tag_data,
118 tag_size);
119 CHKERR mField.get_moab().tag_get_by_ptr(
120 thKappa, &ent, 1, (const void **)&kappaPtr, &kappaSize);
121 }
122 MoFEMFunctionReturnHot(0);
123 }
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition Types.hpp:68
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0

◆ iNitailise()

MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::iNitailise ( const FEMethod *  fe_method)
inline

Initialize history variable data.

Create tag on the prism/interface to store damage history variable

Definition at line 70 of file CohesiveInterfaceElement.hpp.

70 {
71 MoFEMFunctionBeginHot;
72
73 double def_damaged = 0;
74 CHKERR mField.get_moab().tag_get_handle(
75 "DAMAGED_PRISM", 1, MB_TYPE_INTEGER, thDamagedPrism,
76 MB_TAG_CREAT | MB_TAG_SPARSE, &def_damaged);
77 const int def_len = 0;
78 CHKERR mField.get_moab().tag_get_handle(
79 "_KAPPA", def_len, MB_TYPE_DOUBLE, thKappa,
80 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN, nullptr);
81 E0 = youngModulus / h;
82 g0 = ft / E0;
83 kappa1 = 2 * Gf / ft;
84 MoFEMFunctionReturnHot(0);
85 }

◆ updateHistory()

virtual MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::updateHistory ( CommonData &  common_data,
const FEMethod *  fe_method 
)
inlinevirtual

Update history variables when converged.

Definition at line 277 of file CohesiveInterfaceElement.hpp.

278 {
279 MoFEMFunctionBeginHot;
280
281 if (!isInitialised) {
282 CHKERR iNitailise(fe_method);
283 isInitialised = true;
284 }
285 CHKERR getKappa(common_data.gapGlob.size1(), fe_method);
286 bool all_gauss_pts_damaged = true;
287 for (unsigned int gg = 0; gg < common_data.gapGlob.size1(); gg++) {
288 double omega = 0;
289 double g = calcG(gg, common_data.gapLoc);
290 double kappa = fmax(g - g0, kappaPtr[gg]);
291 kappaPtr[gg] = kappa;
292 CHKERR calcOmega(kappa, omega);
293 //if(omega < 1.) {
294 all_gauss_pts_damaged = false;
295 //}
296 }
297 if (all_gauss_pts_damaged) {
298 EntityHandle ent = fe_method->numeredEntFiniteElementPtr->getEnt();
299 int set_prism_as_demaged = 1;
300 CHKERR mField.get_moab().tag_set_data(thDamagedPrism, &ent, 1,
301 &set_prism_as_demaged);
302 }
303 MoFEMFunctionReturnHot(0);
304 }

Member Data Documentation

◆ beta

double CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::beta

Definition at line 59 of file CohesiveInterfaceElement.hpp.

◆ Dglob

MatrixDouble CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::Dglob

Definition at line 125 of file CohesiveInterfaceElement.hpp.

◆ Dloc

MatrixDouble CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::Dloc

Definition at line 125 of file CohesiveInterfaceElement.hpp.

◆ E0

double CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::E0

Definition at line 63 of file CohesiveInterfaceElement.hpp.

◆ ft

double CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::ft

Definition at line 59 of file CohesiveInterfaceElement.hpp.

◆ g0

double CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::g0

Definition at line 63 of file CohesiveInterfaceElement.hpp.

◆ Gf

double CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::Gf

Definition at line 59 of file CohesiveInterfaceElement.hpp.

◆ h

double CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::h

Definition at line 59 of file CohesiveInterfaceElement.hpp.

◆ isInitialised

bool CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::isInitialised

Definition at line 53 of file CohesiveInterfaceElement.hpp.

◆ kappa1

double CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::kappa1

Definition at line 63 of file CohesiveInterfaceElement.hpp.

◆ kappaPtr

double* CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::kappaPtr

Definition at line 99 of file CohesiveInterfaceElement.hpp.

◆ kappaSize

int CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::kappaSize

Definition at line 100 of file CohesiveInterfaceElement.hpp.

◆ mField

MoFEM::Interface& CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::mField

Definition at line 52 of file CohesiveInterfaceElement.hpp.

◆ pRisms

Range CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::pRisms

Definition at line 60 of file CohesiveInterfaceElement.hpp.

◆ thDamagedPrism

Tag CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::thDamagedPrism

Definition at line 61 of file CohesiveInterfaceElement.hpp.

◆ thKappa

Tag CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::thKappa

Definition at line 61 of file CohesiveInterfaceElement.hpp.

◆ youngModulus

double CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::youngModulus

Definition at line 59 of file CohesiveInterfaceElement.hpp.

The documentation for this struct was generated from the following file:
Generated by Doxygen 1.9.8 and hosted at University of Glasgow