v0.14.0
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 
8 namespace CohesiveElement {
9 
10 /** \brief Cohesive element implementation
11 
12  \bug Interface element not working with HO geometry.
13 */
15 
16  struct CommonData {
19  ublas::vector<MatrixDouble > R;
20  };
22 
25  int getRule(int order) { return 2*order; };
26  };
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 
58  virtual ~PhysicalEquation() {
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);
82  E0 = youngModulus/h;
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) {
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 
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  */
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  */
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  */
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  */
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  */
308 
310  FlatPrismElementForcesAndSourcesCore::UserDataOperator(field_name,ForcesAndSourcesCore::UserDataOperator::OPROW) {}
311 
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  */
343 
345  OpCalculateGapGlobal(const std::string field_name,CommonData &common_data):
346  FlatPrismElementForcesAndSourcesCore::UserDataOperator(field_name,ForcesAndSourcesCore::UserDataOperator::OPROW),
347  commonData(common_data) {}
348 
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  */
401 
403  OpCalculateGapLocal(const std::string field_name,CommonData &common_data):
404  FlatPrismElementForcesAndSourcesCore::UserDataOperator(field_name,ForcesAndSourcesCore::UserDataOperator::OPROW),
405  commonData(common_data) {}
406 
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  */
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 
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  */
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 
485  int row_side,int col_side,
486  EntityType row_type,EntityType col_type,
487  EntitiesFieldData::EntData &row_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  */
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 
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));
589  feHistory.getOpPtrVector().push_back(new OpHistory(field_name,commonData,*pit));
590  }
591 
593  }
594 
595 };
596 
597 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::h
double h
Definition: CohesiveInterfaceElement.hpp:61
g
constexpr double g
Definition: shallow_wave.cpp:63
omega
constexpr double omega
Save field DOFS on vertices/tags.
Definition: dynamic_first_order_con_law.cpp:93
CohesiveElement::CohesiveInterfaceElement::OpCalculateGapGlobal::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: CohesiveInterfaceElement.hpp:349
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation
Constitutive (physical) equation for interface.
Definition: CohesiveInterfaceElement.hpp:51
EntityHandle
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::iNitailise
MoFEMErrorCode iNitailise(const FEMethod *fe_method)
Initialize history variable data.
Definition: CohesiveInterfaceElement.hpp:72
CohesiveElement::CohesiveInterfaceElement::feLhs
MyPrism feLhs
Definition: CohesiveInterfaceElement.hpp:28
CohesiveElement::CohesiveInterfaceElement::OpSetSignToShapeFunctions
Set negative sign to shape functions on face 4.
Definition: CohesiveInterfaceElement.hpp:307
CohesiveElement::CohesiveInterfaceElement::OpRhs::physicalEqations
PhysicalEquation & physicalEqations
Definition: CohesiveInterfaceElement.hpp:434
MoFEM::FlatPrismElementForcesAndSourcesCore::FlatPrismElementForcesAndSourcesCore
FlatPrismElementForcesAndSourcesCore(Interface &m_field)
Definition: FlatPrismElementForcesAndSourcesCore.cpp:11
CohesiveElement::CohesiveInterfaceElement::OpLhs::physicalEqations
PhysicalEquation & physicalEqations
Definition: CohesiveInterfaceElement.hpp:478
CohesiveElement::CohesiveInterfaceElement::OpCalculateGapGlobal::OpCalculateGapGlobal
OpCalculateGapGlobal(const std::string field_name, CommonData &common_data)
Definition: CohesiveInterfaceElement.hpp:345
CohesiveElement::CohesiveInterfaceElement::OpRhs
Operator calculate right hand side vector.
Definition: CohesiveInterfaceElement.hpp:431
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
E
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcG
double calcG(int gg, MatrixDouble gap_loc)
Calculate gap opening.
Definition: CohesiveInterfaceElement.hpp:95
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
CohesiveElement::CohesiveInterfaceElement::CommonData
Definition: CohesiveInterfaceElement.hpp:16
CohesiveElement::CohesiveInterfaceElement::OpHistory::physicalEqations
PhysicalEquation & physicalEqations
Definition: CohesiveInterfaceElement.hpp:548
CohesiveElement::CohesiveInterfaceElement::feRhs
MyPrism feRhs
Definition: CohesiveInterfaceElement.hpp:27
CohesiveElement::CohesiveInterfaceElement::OpLhs::K
MatrixDouble K
Definition: CohesiveInterfaceElement.hpp:483
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
CohesiveElement::CohesiveInterfaceElement::commonData
CommonData commonData
Definition: CohesiveInterfaceElement.hpp:21
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
CohesiveElement::CohesiveInterfaceElement::CommonData::gapLoc
MatrixDouble gapLoc
Definition: CohesiveInterfaceElement.hpp:18
CohesiveElement::CohesiveInterfaceElement::OpLhs::D
MatrixDouble D
Definition: CohesiveInterfaceElement.hpp:483
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
order
constexpr int order
Definition: dg_projection.cpp:18
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcTangetDglob
MoFEMErrorCode calcTangetDglob(const double omega, double g, const VectorDouble &gap_loc, MatrixDouble &R)
Calculate tangent material stiffness.
Definition: CohesiveInterfaceElement.hpp:173
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
CohesiveElement::CohesiveInterfaceElement::OpCalculateGapLocal::OpCalculateGapLocal
OpCalculateGapLocal(const std::string field_name, CommonData &common_data)
Definition: CohesiveInterfaceElement.hpp:403
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::thDamagedPrism
Tag thDamagedPrism
Definition: CohesiveInterfaceElement.hpp:63
CohesiveElement::CohesiveInterfaceElement::OpLhs
Operator calculate element stiffens matrix.
Definition: CohesiveInterfaceElement.hpp:475
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
CohesiveElement::CohesiveInterfaceElement::OpCalculateGapLocal::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: CohesiveInterfaceElement.hpp:407
CohesiveElement::CohesiveInterfaceElement::feHistory
MyPrism feHistory
Definition: CohesiveInterfaceElement.hpp:29
FEMethod
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::g0
double g0
Definition: CohesiveInterfaceElement.hpp:65
CohesiveElement::CohesiveInterfaceElement::OpSetSignToShapeFunctions::OpSetSignToShapeFunctions
OpSetSignToShapeFunctions(const std::string field_name)
Definition: CohesiveInterfaceElement.hpp:309
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::isInitialised
bool isInitialised
Definition: CohesiveInterfaceElement.hpp:54
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
CohesiveElement::CohesiveInterfaceElement::getFeLhs
MyPrism & getFeLhs()
Definition: CohesiveInterfaceElement.hpp:39
CohesiveElement::CohesiveInterfaceElement::CommonData::gapGlob
MatrixDouble gapGlob
Definition: CohesiveInterfaceElement.hpp:17
a
constexpr double a
Definition: approx_sphere.cpp:30
CohesiveElement::CohesiveInterfaceElement::OpHistory::OpHistory
OpHistory(const std::string field_name, CommonData &common_data, PhysicalEquation &physical_eqations)
Definition: CohesiveInterfaceElement.hpp:549
R
@ R
Definition: free_surface.cpp:394
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
CohesiveElement::CohesiveInterfaceElement::OpLhs::commonData
CommonData & commonData
Definition: CohesiveInterfaceElement.hpp:477
CohesiveElement::CohesiveInterfaceElement::addOps
MoFEMErrorCode addOps(const std::string field_name, boost::ptr_vector< CohesiveInterfaceElement::PhysicalEquation > &interfaces)
Driver function settting all operators needed for interface element.
Definition: CohesiveInterfaceElement.hpp:568
convert.type
type
Definition: convert.py:64
CohesiveElement::CohesiveInterfaceElement::OpHistory::commonData
CommonData & commonData
Definition: CohesiveInterfaceElement.hpp:547
CohesiveElement::CohesiveInterfaceElement
Cohesive element implementation.
Definition: CohesiveInterfaceElement.hpp:14
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::~PhysicalEquation
virtual ~PhysicalEquation()
Definition: CohesiveInterfaceElement.hpp:58
CohesiveElement::CohesiveInterfaceElement::OpRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: CohesiveInterfaceElement.hpp:440
CohesiveElement::CohesiveInterfaceElement::MyPrism
Definition: CohesiveInterfaceElement.hpp:23
CohesiveElement::CohesiveInterfaceElement::OpLhs::OpLhs
OpLhs(const std::string field_name, CommonData &common_data, PhysicalEquation &physical_eqations)
Definition: CohesiveInterfaceElement.hpp:479
CohesiveElement::CohesiveInterfaceElement::~CohesiveInterfaceElement
virtual ~CohesiveInterfaceElement()
Definition: CohesiveInterfaceElement.hpp:34
CohesiveElement::CohesiveInterfaceElement::OpSetSignToShapeFunctions::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: CohesiveInterfaceElement.hpp:312
CohesiveElement::CohesiveInterfaceElement::OpLhs::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: CohesiveInterfaceElement.hpp:484
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::kappaSize
int kappaSize
Definition: CohesiveInterfaceElement.hpp:100
CohesiveElement::CohesiveInterfaceElement::OpCalculateGapLocal
Operator calculate gap in local coordinate system.
Definition: CohesiveInterfaceElement.hpp:400
CohesiveElement::CohesiveInterfaceElement::getFeRhs
MyPrism & getFeRhs()
Definition: CohesiveInterfaceElement.hpp:38
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcDglob
MoFEMErrorCode calcDglob(const double omega, MatrixDouble &R)
Calculate stiffness material matrix.
Definition: CohesiveInterfaceElement.hpp:136
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcOmega
MoFEMErrorCode calcOmega(const double kappa, double &omega)
Calculate damage.
Definition: CohesiveInterfaceElement.hpp:157
CohesiveElement::CohesiveInterfaceElement::CohesiveInterfaceElement
CohesiveInterfaceElement(MoFEM::Interface &m_field)
Definition: CohesiveInterfaceElement.hpp:31
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::thKappa
Tag thKappa
Definition: CohesiveInterfaceElement.hpp:63
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
kappa
double kappa
Definition: free_surface.cpp:183
Range
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MoFEM::FlatPrismElementForcesAndSourcesCore
FlatPrism finite element.
Definition: FlatPrismElementForcesAndSourcesCore.hpp:27
CohesiveElement::CohesiveInterfaceElement::OpHistory::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: CohesiveInterfaceElement.hpp:553
CohesiveElement::CohesiveInterfaceElement::OpRhs::commonData
CommonData & commonData
Definition: CohesiveInterfaceElement.hpp:433
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::kappaPtr
double * kappaPtr
Definition: CohesiveInterfaceElement.hpp:99
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::ft
double ft
Definition: CohesiveInterfaceElement.hpp:61
CohesiveElement::CohesiveInterfaceElement::getFeHistory
MyPrism & getFeHistory()
Definition: CohesiveInterfaceElement.hpp:40
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
CohesiveElement::CohesiveInterfaceElement::OpCalculateGapGlobal
Operator calculate gap, normal vector and rotation matrix.
Definition: CohesiveInterfaceElement.hpp:342
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calculateTraction
virtual MoFEMErrorCode calculateTraction(VectorDouble &traction, int gg, CommonData &common_data, const FEMethod *fe_method)
Calculate tractions.
Definition: CohesiveInterfaceElement.hpp:203
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::Dglob
MatrixDouble Dglob
Definition: CohesiveInterfaceElement.hpp:122
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::pRisms
Range pRisms
Definition: CohesiveInterfaceElement.hpp:62
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::PhysicalEquation
PhysicalEquation(MoFEM::Interface &m_field)
Definition: CohesiveInterfaceElement.hpp:55
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
CohesiveElement::CohesiveInterfaceElement::OpCalculateGapGlobal::commonData
CommonData & commonData
Definition: CohesiveInterfaceElement.hpp:344
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::youngModulus
double youngModulus
Definition: CohesiveInterfaceElement.hpp:61
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
CohesiveElement::CohesiveInterfaceElement::OpRhs::Nf
VectorDouble Nf
Definition: CohesiveInterfaceElement.hpp:439
CohesiveElement
Definition: arc_length_interface.cpp:29
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
CohesiveElement::CohesiveInterfaceElement::OpLhs::ND
MatrixDouble ND
Definition: CohesiveInterfaceElement.hpp:483
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::updateHistory
virtual MoFEMErrorCode updateHistory(CommonData &common_data, const FEMethod *fe_method)
Update history variables when converged.
Definition: CohesiveInterfaceElement.hpp:272
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
CohesiveElement::CohesiveInterfaceElement::OpCalculateGapLocal::commonData
CommonData & commonData
Definition: CohesiveInterfaceElement.hpp:402
CohesiveElement::CohesiveInterfaceElement::MyPrism::MyPrism
MyPrism(MoFEM::Interface &m_field)
Definition: CohesiveInterfaceElement.hpp:24
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calculateTangentStiffeness
virtual MoFEMErrorCode calculateTangentStiffeness(MatrixDouble &tangent_matrix, int gg, CommonData &common_data, const FEMethod *fe_method)
Calculate tangent stiffness.
Definition: CohesiveInterfaceElement.hpp:231
CohesiveElement::CohesiveInterfaceElement::OpRhs::OpRhs
OpRhs(const std::string field_name, CommonData &common_data, PhysicalEquation &physical_eqations)
Definition: CohesiveInterfaceElement.hpp:435
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::Dloc
MatrixDouble Dloc
Definition: CohesiveInterfaceElement.hpp:122
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::beta
double beta
Definition: CohesiveInterfaceElement.hpp:61
CohesiveElement::CohesiveInterfaceElement::OpRhs::traction
VectorDouble traction
Definition: CohesiveInterfaceElement.hpp:439
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::E0
double E0
Definition: CohesiveInterfaceElement.hpp:65
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
CohesiveElement::CohesiveInterfaceElement::OpHistory
Operator update history variables.
Definition: CohesiveInterfaceElement.hpp:545
CohesiveElement::CohesiveInterfaceElement::MyPrism::getRule
int getRule(int order)
Definition: CohesiveInterfaceElement.hpp:25
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::mField
MoFEM::Interface & mField
Definition: CohesiveInterfaceElement.hpp:53
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::kappa1
double kappa1
Definition: CohesiveInterfaceElement.hpp:65
CohesiveElement::CohesiveInterfaceElement::CommonData::R
ublas::vector< MatrixDouble > R
Definition: CohesiveInterfaceElement.hpp:19
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::Gf
double Gf
Definition: CohesiveInterfaceElement.hpp:61