v0.14.0
Loading...
Searching...
No Matches
CohesiveInterfaceElement.hpp
Go to the documentation of this file.
1/** \file CohesiveInterfaceElement.hpp
2 \brief Implementation of linear interface element
3
4*/
5
6
7
8namespace CohesiveElement {
9
10/** \brief Cohesive element implementation
11
12 \bug Interface element not working with HO geometry.
13*/
15
16 struct CommonData {
17 MatrixDouble gapGlob;
18 MatrixDouble gapLoc;
19 ublas::vector<MatrixDouble > R;
20 };
22
30
32 feRhs(m_field),feLhs(m_field),feHistory(m_field) {};
33
35
36 }
37
38 MyPrism& getFeRhs() { return feRhs; }
39 MyPrism& getFeLhs() { return feLhs; }
41
42 /** \brief Constitutive (physical) equation for interface
43
44 This is linear degradation model. Material parameters are: strength
45 \f$f_t\f$, interface fracture energy \f$G_f\f$, elastic material stiffness
46 \f$E\f$. Parameter \f$\beta\f$ controls how interface opening is calculated.
47
48 Model parameter is interface penalty thickness \f$h\f$.
49
50 */
52
56 mField(m_field),isInitialised(false) {};
57
59 }
60
64
65 double E0,g0,kappa1;
66
67 /** \brief Initialize history variable data
68
69 Create tag on the prism/interface to store damage history variable
70
71 */
72 MoFEMErrorCode iNitailise(const FEMethod *fe_method) {
74
75 double def_damaged = 0;
76 rval = mField.get_moab().tag_get_handle(
77 "DAMAGED_PRISM",1,MB_TYPE_INTEGER,thDamagedPrism,MB_TAG_CREAT|MB_TAG_SPARSE,&def_damaged
78 ); MOAB_THROW(rval);
79 const int def_len = 0;
80 rval = mField.get_moab().tag_get_handle("_KAPPA",def_len,MB_TYPE_DOUBLE,
81 thKappa,MB_TAG_CREAT|MB_TAG_SPARSE|MB_TAG_VARLEN,NULL); CHKERRG(rval);
83 g0 = ft/E0;
84 kappa1 = 2*Gf/ft;
86 }
87
88 /** \brief Calculate gap opening
89
90 \f[
91 g = \sqrt{ g_n^2 + \beta(g_{s1}^2 + g_{s2}^2)}
92 \f]
93
94 */
95 double calcG(int gg,MatrixDouble gap_loc) {
96 return sqrt(pow(gap_loc(gg,0),2)+beta*(pow(gap_loc(gg,1),2)+pow(gap_loc(gg,2),2)));
97 }
98
99 double *kappaPtr;
101
102 /** \brief Get pointer from the mesh to histoy variables \f$\kappa\f$
103 */
104 MoFEMErrorCode getKappa(int nb_gauss_pts,const FEMethod *fe_method) {
106 EntityHandle ent = fe_method->numeredEntFiniteElementPtr->getEnt();
107
108 rval = mField.get_moab().tag_get_by_ptr(thKappa,&ent,1,(const void **)&kappaPtr,&kappaSize);
109 if(rval != MB_SUCCESS || kappaSize != nb_gauss_pts) {
110 VectorDouble kappa;
111 kappa.resize(nb_gauss_pts);
112 kappa.clear();
113 int tag_size[1];
114 tag_size[0] = nb_gauss_pts;
115 void const* tag_data[] = { &kappa[0] };
116 rval = mField.get_moab().tag_set_by_ptr(thKappa,&ent,1,tag_data,tag_size); CHKERRG(rval);
117 rval = mField.get_moab().tag_get_by_ptr(thKappa,&ent,1,(const void **)&kappaPtr,&kappaSize); CHKERRG(rval);
118 }
120 }
121
122 MatrixDouble Dglob,Dloc;
123
124 /** \brief Calculate stiffness material matrix
125
126 \f[
127 \mathbf{D}_\textrm{loc} = (1-\Omega) \mathbf{I} E_0
128 \f]
129 where \f$E_0\f$ is initial interface penalty stiffness
130
131 \f[
132 \mathbf{D}_\textrm{glob} = \mathbf{R}^\textrm{T} \mathbf{D}_\textrm{loc}\mathbf{R}
133 \f]
134
135 */
136 MoFEMErrorCode calcDglob(const double omega,MatrixDouble &R) {
138 Dglob.resize(3,3);
139 Dloc.resize(3,3);
140 Dloc.clear();
141 double E = (1-omega)*E0;
142 Dloc(0,0) = E;
143 Dloc(1,1) = E;
144 Dloc(2,2) = E;
145 Dglob = prod( Dloc, R );
146 Dglob = prod( trans(R), Dglob );
148 }
149
150 /** \brief Calculate damage
151
152 \f[
153 \Omega = \frac{1}{2} \frac{(2 G_f E_0+f_t^2)\kappa}{(ft+E_0 \kappa)G_f}
154 \f]
155
156 */
157 MoFEMErrorCode calcOmega(const double kappa,double& omega) {
159 omega = 0;
160 if(kappa>=kappa1) {
161 omega = 1;
163 } else if(kappa>0) {
164 double a = (2.0*Gf*E0+ft*ft)*kappa;
165 double b = (ft+E0*kappa)*Gf;
166 omega = 0.5*a/b;
167 }
169 }
170
171 /** \brief Calculate tangent material stiffness
172 */
173 MoFEMErrorCode calcTangetDglob(const double omega,double g,const VectorDouble& gap_loc,MatrixDouble &R) {
175 Dglob.resize(3,3);
176 Dloc.resize(3,3);
177 double domega = 0.5*(2*Gf*E0+ft*ft)/((ft+(g-ft/E0)*E0)*Gf) - 0.5*((g-ft/E0)*(2*Gf*E0+ft*ft)*E0)/(pow(ft+(g-ft/E0)*E0,2)*Gf);
178 Dloc.resize(3,3);
179 //r0
180 Dloc(0,0) = (1-omega)*E0 - domega*E0*gap_loc[0]*gap_loc[0]/g;
181 Dloc(0,1) = -domega*E0*gap_loc[0]*beta*gap_loc[1]/g;
182 Dloc(0,2) = -domega*E0*gap_loc[0]*beta*gap_loc[2]/g;
183 //r1
184 Dloc(1,0) = -domega*E0*gap_loc[1]*gap_loc[0]/g;
185 Dloc(1,1) = (1-omega)*E0 - domega*E0*gap_loc[1]*beta*gap_loc[1]/g;
186 Dloc(1,2) = -domega*E0*gap_loc[1]*beta*gap_loc[2]/g;
187 //r2
188 Dloc(2,0) = -domega*E0*gap_loc[2]*gap_loc[0]/g;
189 Dloc(2,1) = -domega*E0*gap_loc[2]*beta*gap_loc[1]/g;
190 Dloc(2,2) = (1-omega)*E0 - domega*E0*gap_loc[2]*beta*gap_loc[2]/g;
191 Dglob = prod(Dloc,R);
192 Dglob = prod(trans(R),Dglob);
194 }
195
196 /** \brief Calculate tractions
197
198 \f[
199 \mathbf{t} = \mathbf{D}_\textrm{glob}\mathbf{g}
200 \f]
201
202 */
203 virtual MoFEMErrorCode calculateTraction(
204 VectorDouble &traction,
205 int gg,CommonData &common_data,
206 const FEMethod *fe_method
207 ) {
209
210 if(!isInitialised) {
211 ierr = iNitailise(fe_method); CHKERRG(ierr);
212 isInitialised = true;
213 }
214 if(gg==0) {
215 ierr = getKappa(common_data.gapGlob.size1(),fe_method); CHKERRG(ierr);
216 }
217 double g = calcG(gg,common_data.gapLoc);
218 double kappa = fmax(g-g0,kappaPtr[gg]);
219 double omega = 0;
221 //std::cerr << gg << " " << omega << std::endl;
222 ierr = calcDglob(omega,common_data.R[gg]); CHKERRG(ierr);
223 traction.resize(3);
224 ublas::matrix_row<MatrixDouble > gap_glob(common_data.gapGlob,gg);
225 noalias(traction) = prod(Dglob,gap_glob);
227 }
228
229 /** \brief Calculate tangent stiffness
230 */
231 virtual MoFEMErrorCode calculateTangentStiffeness(
232 MatrixDouble &tangent_matrix,
233 int gg,CommonData &common_data,
234 const FEMethod *fe_method
235 ) {
237
238 try {
239 if(!isInitialised) {
240 ierr = iNitailise(fe_method); CHKERRG(ierr);
241 isInitialised = true;
242 }
243 if(gg==0) {
244 ierr = getKappa(common_data.gapGlob.size1(),fe_method); CHKERRG(ierr);
245 }
246 double g = calcG(gg,common_data.gapLoc);
247 double kappa = fmax(g-g0,kappaPtr[gg]);
248 double omega = 0;
250 //std::cerr << gg << " " << omega << std::endl;
251 int iter;
252 ierr = SNESGetIterationNumber(fe_method->snes,&iter); CHKERRG(ierr);
253 if((kappa <= kappaPtr[gg])||(kappa>=kappa1)||(iter <= 1)) {
254 ierr = calcDglob(omega,common_data.R[gg]); CHKERRG(ierr);
255 } else {
256 ublas::matrix_row<MatrixDouble > g_loc(common_data.gapLoc,gg);
257 ierr = calcTangetDglob(omega,g,g_loc,common_data.R[gg]);
258 }
259 tangent_matrix.resize(3,3);
260 noalias(tangent_matrix) = Dglob;
261 //std::cerr << "t " << tangent_matrix << std::endl;
262 } catch (const std::exception& ex) {
263 std::ostringstream ss;
264 ss << "throw in method: " << ex.what() << std::endl;
265 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
266 }
268 }
269
270 /** \brief Update history variables when converged
271 */
272 virtual MoFEMErrorCode updateHistory(
273 CommonData &common_data,const FEMethod *fe_method
274 ) {
276
277
278 if(!isInitialised) {
279 ierr = iNitailise(fe_method); CHKERRG(ierr);
280 isInitialised = true;
281 }
282 ierr = getKappa(common_data.gapGlob.size1(),fe_method); CHKERRG(ierr);
283 bool all_gauss_pts_damaged = true;
284 for(unsigned int gg = 0;gg<common_data.gapGlob.size1();gg++) {
285 double omega = 0;
286 double g = calcG(gg,common_data.gapLoc);
287 double kappa = fmax(g-g0,kappaPtr[gg]);
288 kappaPtr[gg] = kappa;
290 //if(omega < 1.) {
291 all_gauss_pts_damaged = false;
292 //}
293 }
294 if(all_gauss_pts_damaged) {
295 EntityHandle ent = fe_method->numeredEntFiniteElementPtr->getEnt();
296 int set_prism_as_demaged = 1;
297 rval = mField.get_moab().tag_set_data(thDamagedPrism,&ent,1,&set_prism_as_demaged); CHKERRG(rval);
298 }
300 }
301
302 };
303
304
305 /** \brief Set negative sign to shape functions on face 4
306 */
307 struct OpSetSignToShapeFunctions: public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
308
310 FlatPrismElementForcesAndSourcesCore::UserDataOperator(field_name,ForcesAndSourcesCore::UserDataOperator::OPROW) {}
311
312 MoFEMErrorCode doWork(int side,EntityType type,EntitiesFieldData::EntData &data) {
314 if(data.getN().size1()==0) MoFEMFunctionReturnHot(0);
315 if(data.getN().size2()==0) MoFEMFunctionReturnHot(0);
316 switch(type) {
317 case MBVERTEX:
318 for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
319 for(int nn = 3;nn<6;nn++) {
320 data.getN()(gg,nn) *= -1;
321 }
322 }
323 break;
324 case MBEDGE:
325 if(side < 3) MoFEMFunctionReturnHot(0);
326 data.getN() *= -1;
327 break;
328 case MBTRI:
329 if(side == 3) MoFEMFunctionReturnHot(0);
330 data.getN() *= -1;
331 break;
332 default:
333 SETERRQ(PETSC_COMM_SELF,1,"data inconsitency");
334 }
336 }
337
338 };
339
340 /** \brief Operator calculate gap, normal vector and rotation matrix
341 */
342 struct OpCalculateGapGlobal: public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
343
345 OpCalculateGapGlobal(const std::string field_name,CommonData &common_data):
346 FlatPrismElementForcesAndSourcesCore::UserDataOperator(field_name,ForcesAndSourcesCore::UserDataOperator::OPROW),
347 commonData(common_data) {}
348
349 MoFEMErrorCode doWork(
350 int side,EntityType type,EntitiesFieldData::EntData &data) {
352 try {
353 int nb_dofs = data.getIndices().size();
354 if(nb_dofs == 0) MoFEMFunctionReturnHot(0);
355 int nb_gauss_pts = data.getN().size1();
356 if(type == MBVERTEX) {
357 commonData.R.resize(nb_gauss_pts);
358 for(int gg = 0;gg<nb_gauss_pts;gg++) {
359 commonData.R[gg].resize(3,3);
360 double nrm2_normal = 0;
361 double nrm2_tangent1 = 0;
362 double nrm2_tangent2 = 0;
363 for(int dd = 0;dd<3;dd++) {
364 nrm2_normal += pow(getNormalsAtGaussPtsF3()(gg,dd),2);
365 nrm2_tangent1 += pow(getTangent1AtGaussPtF3()(gg,dd),2);
366 nrm2_tangent2 += pow(getTangent2AtGaussPtF3()(gg,dd),2);
367 }
368 nrm2_normal = sqrt(nrm2_normal);
369 nrm2_tangent1 = sqrt(nrm2_tangent1);
370 nrm2_tangent2 = sqrt(nrm2_tangent2);
371 for(int dd = 0;dd<3;dd++) {
372 commonData.R[gg](0,dd) = getNormalsAtGaussPtsF3()(gg,dd)/nrm2_normal;
373 commonData.R[gg](1,dd) = getTangent1AtGaussPtF3()(gg,dd)/nrm2_tangent1;
374 commonData.R[gg](2,dd) = getTangent2AtGaussPtF3()(gg,dd)/nrm2_tangent2;
375 }
376 }
377 }
378 if(type == MBVERTEX) {
379 commonData.gapGlob.resize(nb_gauss_pts,3);
380 commonData.gapGlob.clear();
381 }
382 for(int gg = 0;gg<nb_gauss_pts;gg++) {
383 for(int dd = 0;dd<3;dd++) {
384 commonData.gapGlob(gg,dd) += cblas_ddot(
385 nb_dofs/3,&data.getN(gg)[0],1,&data.getFieldData()[dd],3);
386 }
387 }
388 } catch (const std::exception& ex) {
389 std::ostringstream ss;
390 ss << "throw in method: " << ex.what() << std::endl;
391 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
392 }
394 }
395
396 };
397
398 /** \brief Operator calculate gap in local coordinate system
399 */
400 struct OpCalculateGapLocal: public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
401
403 OpCalculateGapLocal(const std::string field_name,CommonData &common_data):
404 FlatPrismElementForcesAndSourcesCore::UserDataOperator(field_name,ForcesAndSourcesCore::UserDataOperator::OPROW),
405 commonData(common_data) {}
406
407 MoFEMErrorCode doWork(int side,EntityType type,EntitiesFieldData::EntData &data) {
409 try {
410 if(type == MBVERTEX) {
411 int nb_gauss_pts = data.getN().size1();
412 commonData.gapLoc.resize(nb_gauss_pts,3);
413 for(int gg = 0;gg<nb_gauss_pts;gg++) {
414 ublas::matrix_row<MatrixDouble > gap_glob(commonData.gapGlob,gg);
415 ublas::matrix_row<MatrixDouble > gap_loc(commonData.gapLoc,gg);
416 gap_loc = prod(commonData.R[gg],gap_glob);
417 }
418 }
419 } catch (const std::exception& ex) {
420 std::ostringstream ss;
421 ss << "throw in method: " << ex.what() << std::endl;
422 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
423 }
425 }
426
427 };
428
429 /** \brief Operator calculate right hand side vector
430 */
431 struct OpRhs: public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
432
435 OpRhs(const std::string field_name,CommonData &common_data,PhysicalEquation &physical_eqations):
436 FlatPrismElementForcesAndSourcesCore::UserDataOperator(field_name,ForcesAndSourcesCore::UserDataOperator::OPROW),
437 commonData(common_data),physicalEqations(physical_eqations) {}
438
439 VectorDouble traction,Nf;
440 MoFEMErrorCode doWork(int side,EntityType type,EntitiesFieldData::EntData &data) {
442
443 try {
444 int nb_dofs = data.getIndices().size();
445 if(nb_dofs == 0) MoFEMFunctionReturnHot(0);
446 if(physicalEqations.pRisms.find(getNumeredEntFiniteElementPtr()->getEnt()) == physicalEqations.pRisms.end()) {
448 }
449 Nf.resize(nb_dofs);
450 Nf.clear();
451 int nb_gauss_pts = data.getN().size1();
452 for(int gg = 0;gg<nb_gauss_pts;gg++) {
454 double w = getGaussPts()(2,gg)*cblas_dnrm2(3,&getNormalsAtGaussPtsF3()(gg,0),1)*0.5;
455 for(int nn = 0;nn<nb_dofs/3;nn++) {
456 for(int dd = 0;dd<3;dd++) {
457 Nf[3*nn+dd] += w*data.getN(gg)[nn]*traction[dd];
458 }
459 }
460 }
461 ierr = VecSetValues(getFEMethod()->snes_f,
462 data.getIndices().size(),&data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRG(ierr);
463 } catch (const std::exception& ex) {
464 std::ostringstream ss;
465 ss << "throw in method: " << ex.what() << std::endl;
466 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
467 }
469 }
470
471 };
472
473 /** \brief Operator calculate element stiffens matrix
474 */
475 struct OpLhs: public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
476
479 OpLhs(const std::string field_name,CommonData &common_data,PhysicalEquation &physical_eqations):
480 FlatPrismElementForcesAndSourcesCore::UserDataOperator(field_name,ForcesAndSourcesCore::UserDataOperator::OPROWCOL),
481 commonData(common_data),physicalEqations(physical_eqations) { sYmm = false; }
482
483 MatrixDouble K,D,ND;
484 MoFEMErrorCode doWork(
485 int row_side,int col_side,
486 EntityType row_type,EntityType col_type,
487 EntitiesFieldData::EntData &row_data,
488 EntitiesFieldData::EntData &col_data
489 ) {
491
492 try {
493 int nb_row = row_data.getIndices().size();
494 if(nb_row == 0) MoFEMFunctionReturnHot(0);
495 int nb_col = col_data.getIndices().size();
496 if(nb_col == 0) MoFEMFunctionReturnHot(0);
497 if(physicalEqations.pRisms.find(getNumeredEntFiniteElementPtr()->getEnt())
498 == physicalEqations.pRisms.end()) {
500 }
501 //std::cerr << row_side << " " << row_type << " " << row_data.getN() << std::endl;
502 //std::cerr << col_side << " " << col_type << " " << col_data.getN() << std::endl;
503 ND.resize(nb_row,3);
504 K.resize(nb_row,nb_col);
505 K.clear();
506 int nb_gauss_pts = row_data.getN().size1();
507 for(int gg = 0;gg<nb_gauss_pts;gg++) {
509 double w = getGaussPts()(2,gg)*cblas_dnrm2(3,&getNormalsAtGaussPtsF3()(gg,0),1)*0.5;
510 ND.clear();
511 for(int nn = 0; nn<nb_row/3;nn++) {
512 for(int dd = 0;dd<3;dd++) {
513 for(int DD = 0;DD<3;DD++) {
514 ND(3*nn+dd,DD) += row_data.getN(gg)[nn]*D(dd,DD);
515 }
516 }
517 }
518 for(int nn = 0; nn<nb_row/3; nn++) {
519 for(int dd = 0;dd<3;dd++) {
520 for(int NN = 0; NN<nb_col/3; NN++) {
521 for(int DD = 0; DD<3;DD++) {
522 K(3*nn+dd,3*NN+DD) += w*ND(3*nn+dd,DD)*col_data.getN(gg)[NN];
523 }
524 }
525 }
526 }
527 }
528 ierr = MatSetValues(getFEMethod()->snes_B,
529 nb_row,&row_data.getIndices()[0],
530 nb_col,&col_data.getIndices()[0],
531 &K(0,0),ADD_VALUES
532 ); CHKERRG(ierr);
533 } catch (const std::exception& ex) {
534 std::ostringstream ss;
535 ss << "throw in method: " << ex.what() << std::endl;
536 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
537 }
539 }
540
541 };
542
543 /** \brief Operator update history variables
544 */
545 struct OpHistory: public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
546
549 OpHistory(const std::string field_name,CommonData &common_data,PhysicalEquation &physical_eqations):
550 FlatPrismElementForcesAndSourcesCore::UserDataOperator(field_name,ForcesAndSourcesCore::UserDataOperator::OPROW),
551 commonData(common_data),physicalEqations(physical_eqations) {}
552
553 MoFEMErrorCode doWork(int side,EntityType type,EntitiesFieldData::EntData &data) {
555
556 if(type != MBVERTEX) MoFEMFunctionReturnHot(0);
557 if(physicalEqations.pRisms.find(getNumeredEntFiniteElementPtr()->getEnt()) == physicalEqations.pRisms.end()) {
559 }
562 }
563
564 };
565
566 /** \brief Driver function settting all operators needed for interface element
567 */
568 MoFEMErrorCode addOps(const std::string field_name,boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation> &interfaces) {
570
571 //Rhs
575 //Lhs
579 //History
583
584 //add equations/data for physical interfaces
585 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator pit;
586 for(pit = interfaces.begin();pit!=interfaces.end();pit++) {
587 feRhs.getOpPtrVector().push_back(new OpRhs(field_name,commonData,*pit));
588 feLhs.getOpPtrVector().push_back(new OpLhs(field_name,commonData,*pit));
590 }
591
593 }
594
595};
596
597}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
static PetscErrorCode ierr
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
constexpr double omega
Save field DOFS on vertices/tags.
double kappa
@ R
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
constexpr auto field_name
constexpr double g
Operator calculate gap, normal vector and rotation matrix.
OpCalculateGapGlobal(const std::string field_name, CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpCalculateGapLocal(const std::string field_name, CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpHistory(const std::string field_name, CommonData &common_data, PhysicalEquation &physical_eqations)
Operator calculate element stiffens matrix.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpLhs(const std::string field_name, CommonData &common_data, PhysicalEquation &physical_eqations)
OpRhs(const std::string field_name, CommonData &common_data, PhysicalEquation &physical_eqations)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
double calcG(int gg, MatrixDouble gap_loc)
Calculate gap opening.
virtual MoFEMErrorCode calculateTangentStiffeness(MatrixDouble &tangent_matrix, int gg, CommonData &common_data, const FEMethod *fe_method)
Calculate tangent stiffness.
MoFEMErrorCode iNitailise(const FEMethod *fe_method)
Initialize history variable data.
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.
MoFEMErrorCode calcDglob(const double omega, MatrixDouble &R)
Calculate stiffness material matrix.
virtual MoFEMErrorCode updateHistory(CommonData &common_data, const FEMethod *fe_method)
Update history variables when converged.
MoFEMErrorCode getKappa(int nb_gauss_pts, const FEMethod *fe_method)
Get pointer from the mesh to histoy variables .
virtual MoFEMErrorCode calculateTraction(VectorDouble &traction, int gg, CommonData &common_data, const FEMethod *fe_method)
Calculate tractions.
MoFEMErrorCode addOps(const std::string field_name, boost::ptr_vector< CohesiveInterfaceElement::PhysicalEquation > &interfaces)
Driver function settting all operators needed for interface element.
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.