v0.14.0
Loading...
Searching...
No Matches
NitschePeriodicMethod.hpp
Go to the documentation of this file.
1/** \file NitschePeriodicMethod.hpp
2 * \ingroup nitsche_method
3 * \brief Implementation of Nitsche's method for periodic constrains.
4 */
5
6/*
7 * This file is part of MoFEM.
8 * MoFEM is free software: you can redistribute it and/or modify it under
9 * the terms of the GNU Lesser General Public License as published by the
10 * Free Software Foundation, either version 3 of the License, or (at your
11 * option) any later version.
12 *
13 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20
21#ifndef __NITSCHE_PERIODIC_METHOD_HPP__
22#define __NITSCHE_PERIODIC_METHOD_HPP__
23
24/** \brief Periodic Nitsche method
25* \ingroup nitsche_method
26*/
28
29 struct CommonData {
30
33 double eps;
34
35 map<EntityHandle,vector<VectorDouble > > localCoordsMap;
36 map<EntityHandle,vector<int> > inTetFaceGaussPtsNumber;
37 map<EntityHandle,vector<int> > inTetTetGaussPtsNumber;
38 map<EntityHandle,double> dIstance;
39 vector<VectorDouble > dIsplacements;
40 vector<VectorDouble > coordsAtGaussPts;
41 vector<VectorDouble > hoCoordsAtGaussPts;
42 vector<MatrixDouble > stressJacobian;
43 vector<MatrixDouble > sTress;
44
46 int gG;
47 int sIde;
49 ublas::vector<int> iNdices;
50 ublas::vector<int> dofOrders;
51 VectorDouble shapeFunctions;
52 MatrixDouble diffShapeFunctions;
54 int gg,int side,EntityType type
55 ):
56 gG(gg),
57 sIde(side),
58 tYpe(type) {
59 }
60 };
61
62 typedef multi_index_container<
63 MultiIndexData,
64 indexed_by<
65 ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData,int,MultiIndexData::gG) >,
66 ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData,int,MultiIndexData::sIde) >,
67 ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData,EntityType,MultiIndexData::tYpe) >,
68 ordered_unique<
69 composite_key<
70 MultiIndexData,
71 member<MultiIndexData,int,&MultiIndexData::gG>,
72 member<MultiIndexData,int,&MultiIndexData::sIde>,
73 member<MultiIndexData,EntityType,&MultiIndexData::tYpe>
74 > >
75 >
79
80 };
81
82
83 struct OpGetFaceData: FaceElementForcesAndSourcesCore::UserDataOperator {
84
87
88 OpGetFaceData(CommonData &common_data,bool field_disp = true):
90 commonData(common_data),
91 fieldDisp(field_disp) {
92 }
93
94 PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
95 PetscFunctionBegin;
96
97 if(data.getIndices().size()==0) PetscFunctionReturn(0);
98 try {
99 EntityHandle this_face = getNumeredEntFiniteElementPtr()->getEnt();
100 int nb_gauss_pts_on_this_face = data.getN().size1();
101 for(int ffgg = 0;ffgg<nb_gauss_pts_on_this_face;ffgg++) {
102 int gg = commonData.inTetFaceGaussPtsNumber[this_face][ffgg];
103 CommonData::MultiIndexData gauss_pt_data(gg,side,type);
104 pair<CommonData::Container::iterator,bool> p;
105 p = commonData.facesContainer.insert(gauss_pt_data);
106 if(!p.second) {
107 SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"data not inserted");
108 }
109 CommonData::MultiIndexData &p_data = const_cast<CommonData::MultiIndexData&>(*p.first);
110 VectorDouble &shape_fun = p_data.shapeFunctions;
111 int nb_shape_fun = data.getN().size2();
112 shape_fun.resize(nb_shape_fun);
113 cblas_dcopy(nb_shape_fun,&data.getN()(ffgg,0),1,&shape_fun[0],1);
114 p_data.iNdices = data.getIndices();
115 p_data.dofOrders.resize(data.getFieldDofs().size(),false);
116 for(unsigned int dd = 0;dd<data.getFieldDofs().size();dd++) {
117 p_data.dofOrders[dd] = data.getFieldDofs()[dd]->getDofOrder();
118 }
119 int nb_dofs = data.getFieldData().size();
120 if(type == MBVERTEX) {
121 commonData.dIsplacements[gg].resize(3,false);
122 commonData.coordsAtGaussPts[gg].resize(3,false);
123 commonData.hoCoordsAtGaussPts[gg].resize(3,false);
124 for(int rr = 0;rr<3;rr++) {
125 (commonData.dIsplacements[gg])[rr] = cblas_ddot(
126 nb_dofs/3,&data.getFieldData()[rr],3,&data.getN()(ffgg,0),1
127 );
128 if(!fieldDisp) {
129 (commonData.dIsplacements[gg])[rr] -= getCoordsAtGaussPts()(ffgg,rr);
130 }
131 commonData.coordsAtGaussPts[gg][rr] = getCoordsAtGaussPts()(ffgg,rr);
132 commonData.hoCoordsAtGaussPts[gg][rr] = getCoordsAtGaussPts()(ffgg,rr);
133 }
134 } else {
135 for(int rr = 0;rr<3;rr++) {
136 commonData.dIsplacements[gg][rr] += cblas_ddot(
137 nb_dofs/3,&data.getFieldData()[rr],3,&data.getN()(ffgg,0),1
138 );
139 }
140 }
141 }
142 } catch (const std::exception& ex) {
143 ostringstream ss;
144 ss << "throw in method: " << ex.what() << endl;
145 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
146 }
147
148 //cerr << "done\n";
149 PetscFunctionReturn(0);
150 }
151
152 };
153
155
159 commonData(common_data) {}
160
161 int getRule(int order) { return -1; }
162
163 PetscErrorCode setGaussPts(int order) {
164 PetscFunctionBegin;
165
166 try {
167
168 EntityHandle face = numeredEntFiniteElementPtr->getEnt();
169 gaussPts.resize(3,commonData.localCoordsMap[face].size());
170 //cerr << commonData.localCoordsMap[face].size() << endl;
171 for(unsigned int ffgg = 0;ffgg!=commonData.localCoordsMap[face].size();ffgg++) {
172 //cerr << commonData.localCoordsMap[face][ffgg] << endl;
173 gaussPts(0,ffgg) = commonData.localCoordsMap[face][ffgg][0];
174 gaussPts(1,ffgg) = commonData.localCoordsMap[face][ffgg][1];
175 gaussPts(2,ffgg) = 0; // not used, since filed on this face is not integrated
176 }
177
178 } catch (const std::exception& ex) {
179 ostringstream ss;
180 ss << "throw in method: " << ex.what() << endl;
181 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
182 }
183 PetscFunctionReturn(0);
184 }
185
186 };
187
188 struct OpGetVolumeData: VolumeElementForcesAndSourcesCore::UserDataOperator {
189
192
195 CommonData &periodic_common_data
196 ):
197 VolumeElementForcesAndSourcesCore::UserDataOperator("DISPLACEMENT",OPROW),
198 commonData(common_data),
199 periodicCommonData(periodic_common_data) {
200 }
201
202 PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
203 PetscFunctionBegin;
204 if(data.getIndices().size()==0) PetscFunctionReturn(0);
205 try {
206 EntityHandle this_tet = getNumeredEntFiniteElementPtr()->getEnt();
207 int nb_gauss_pts = data.getN().size1();
208 for(int ffgg = 0;ffgg<nb_gauss_pts;ffgg++) {
209 int gg = periodicCommonData.inTetTetGaussPtsNumber[this_tet][ffgg];
210 CommonData::MultiIndexData gauss_pt_data(gg,side,type);
211 pair<CommonData::Container::iterator,bool> p;
212 p = periodicCommonData.volumesContainer.insert(gauss_pt_data);
213 if(!p.second) {
214 SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"data not inserted");
215 }
216 CommonData::MultiIndexData &p_data = const_cast<CommonData::MultiIndexData&>(*p.first);
217 MatrixDouble &diff_shape_fun = p_data.diffShapeFunctions;
218 diff_shape_fun = data.getDiffN(ffgg);
219 p_data.iNdices = data.getIndices();
220 if(type == MBVERTEX) {
223 }
224 }
225 } catch (const std::exception& ex) {
226 ostringstream ss;
227 ss << "throw in method: " << ex.what() << endl;
228 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
229 }
230 PetscFunctionReturn(0);
231 }
232
233 };
234
235 struct PeriodicVolume: VolumeElementForcesAndSourcesCore {
236
239
241 VolumeElementForcesAndSourcesCore(m_field),
242 forFace(0),
243 commonData(common_data) {}
244 int getRule(int order) { return -1; }
245 PetscErrorCode setGaussPts(int order) {
246
247 PetscFunctionBegin;
248 try {
249 gaussPts.resize(4,0,false);
250 EntityHandle tet = numeredEntFiniteElementPtr->getEnt();
251 int vffgg = 0;
252 for(int ff = 0;ff<4;ff++) {
253 EntityHandle face;
254 rval = mField.get_moab().side_element(tet,2,ff,face); CHKERRQ_MOAB(rval);
255 if(face!=forFace) continue;
256 if(commonData.localCoordsMap.find(face)!=commonData.localCoordsMap.end()) {
257 const double coords_tet[12] = { 0,0,0, 1,0,0, 0,1,0, 0,0,1 };
258 int nb_face_gauss_pts = commonData.localCoordsMap[face].size();
259 gaussPts.resize(4,nb_face_gauss_pts);
260 commonData.inTetTetGaussPtsNumber[tet].resize(nb_face_gauss_pts);
261 for(int ffgg = 0;ffgg!=nb_face_gauss_pts;ffgg++,vffgg++) {
263 double N[3];
264 VectorDouble &local_coords = commonData.localCoordsMap[face][ffgg];
265 ierr = ShapeMBTRI(N,&local_coords[0],&local_coords[1],1); CHKERRQ(ierr);
266 for(int dd = 0;dd<3;dd++) {
267 gaussPts(dd,vffgg) =
268 N[0]*coords_tet[3*dataH1.facesNodes(ff,0)+dd]+
269 N[1]*coords_tet[3*dataH1.facesNodes(ff,1)+dd]+
270 N[2]*coords_tet[3*dataH1.facesNodes(ff,2)+dd];
271 }
272 gaussPts(3,vffgg) = 0;
273 }
274 } else {
275 SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"data inconstancy");
276 }
277 }
278 //cerr << gaussPts << endl;
279 } catch (const std::exception& ex) {
280 ostringstream ss;
281 ss << "throw in method: " << ex.what() << endl;
282 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
283 }
284 PetscFunctionReturn(0);
285 }
286
287 };
288
290
294
295 bool iNit;
296 AdaptiveKDTree tRee;
297
299 MoFEM::Interface &m_field,
300 NitscheMethod::BlockData &block_data,
301 NitscheMethod::CommonData &common_data,
302 CommonData &periodic_common_data
303 ):
304 NitscheMethod::MyVolumeFE(m_field,block_data,common_data),
305 periodicCommonData(periodic_common_data),
306 periodicFace(m_field,periodic_common_data),
307 periodicVolume(m_field,periodic_common_data),
308 iNit(false),
309 tRee(&m_field.get_moab()) {
310 periodicFace.getOpPtrVector().push_back(
311 new OpGetFaceData(periodic_common_data)
312 );
313 }
314
316 PetscErrorCode iNitialise() {
317 PetscFunctionBegin;
318
319 treeRootSet = 0;
320 rval = tRee.build_tree(periodicCommonData.skinFaces,&treeRootSet); CHKERRQ_MOAB(rval);
321 //cerr << blockData.fAces << endl;
322 PetscFunctionReturn(0);
323 }
324
326 PetscFunctionBegin;
327
328
329
330
331 if(!iNit) {
332 ierr = iNitialise(); CHKERRQ(ierr);
333 iNit = true;
334 }
335
336 //cerr << "part1\n";
337 try {
338 const double tol = 1e-18;
350 int gg = 0;
351 for(int ff = 0;ff<4;ff++) {
352 if(commonData.facesFePtr[ff]==NULL) continue;
353 for(unsigned int fgg = 0;fgg<commonData.faceGaussPts[ff].size2();fgg++,gg++) {
354 VectorDouble& ray = commonData.rAy[ff];
355 VectorAdaptor coords = VectorAdaptor(
356 3,ublas::shallow_array_adaptor<double>(
357 3,&commonData.hoCoordsAtGaussPts[ff](fgg,0)
358 )
359 );
360 vector<EntityHandle> triangles_out;
361 vector<double> distance_out;
362 rval = tRee.ray_intersect_triangles(
363 treeRootSet,tol,&ray[0],&coords[0],triangles_out,distance_out
364 ); CHKERRQ_MOAB(rval);
365 double diffNTRI[6];
366 ierr = ShapeDiffMBTRI(diffNTRI); CHKERRQ(ierr);
367 vector<EntityHandle> triangles_out_on_other_side;
368 vector<double> distance_out_on_other_side;
369 {
370 vector<EntityHandle>::iterator hit = triangles_out.begin();
371 vector<double>::iterator dit = distance_out.begin();
372 for(;dit!=distance_out.end();dit++,hit++) {
373 const EntityHandle* conn;
374 int number_nodes;
375 rval = mField.get_moab().get_connectivity(*hit,conn,number_nodes,true); CHKERRQ_MOAB(rval);
376 double coords[9];
377 rval = mField.get_moab().get_coords(conn,number_nodes,coords); CHKERRQ_MOAB(rval);
378 double normal[3];
379 ierr = ShapeFaceNormalMBTRI(diffNTRI,coords,normal); CHKERRQ(ierr);
380 *dit *= cblas_ddot(3,normal,1,&ray[0],1)/cblas_dnrm2(3,normal,1);
381 if(*dit>tol) {
382 triangles_out_on_other_side.push_back(*hit);
383 distance_out_on_other_side.push_back(*dit);
384 }
385 }
386 }
387 vector<double>::iterator max_dit = max_element(
388 distance_out_on_other_side.begin(),distance_out_on_other_side.end()
389 );
390 EntityHandle other_face = triangles_out_on_other_side[
391 std::distance(distance_out_on_other_side.begin(),max_dit)
392 ];
393 VectorDouble other_coords = coords+ray*(*max_dit);
394 periodicCommonData.dIstance[other_face] = *max_dit;
395 VectorDouble local_coords;
396 local_coords.resize(2);
397 {
398 const EntityHandle* conn;
399 int number_nodes;
400 rval = mField.get_moab().get_connectivity(other_face,conn,number_nodes,false); CHKERRQ_MOAB(rval);
401
402 {
403 double face_coords[18];
404 rval = mField.get_moab().get_coords(conn,number_nodes,face_coords); CHKERRQ_MOAB(rval);
405
406 VectorDouble res,d_local_coords;
407 res.resize(3,false);
408 d_local_coords.resize(2,false);
409
410 double N[number_nodes];
411 double diffNTRI[number_nodes*2];
412 double nrm2,nrm2_res;
413 const double eps = 1e-14;
414 const int max_it = 10;
415
416 int it = 0;
417 local_coords.clear();
418 do {
419 if(number_nodes == 3) {
420 ierr = ShapeMBTRI(N,&local_coords[0],&local_coords[1],1); CHKERRQ(ierr);
421 } else if(number_nodes == 6) {
422 ierr = ShapeMBTRIQ(N,&local_coords[0],&local_coords[1],1); CHKERRQ(ierr);
423 } else {
424 SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"not implemented");
425 }
426 res.clear();
427 for(int dd = 0;dd<3;dd++) {
428 res[dd] += cblas_ddot(number_nodes,N,1,&face_coords[dd],3)-other_coords[dd];
429 }
430 nrm2_res = cblas_dnrm2(3,&res[0],1);
431 if(number_nodes == 3) {
432 ierr = ShapeDiffMBTRI(diffNTRI); CHKERRQ(ierr);
433 } else {
434 ierr = ShapeDiffMBTRIQ(diffNTRI,&local_coords[0],&local_coords[1],1); CHKERRQ(ierr);
435 }
436 MatrixDouble A;
437 A.resize(3,2);
438 for(int dd = 0;dd<3;dd++) {
439 A(dd,0) = cblas_ddot(number_nodes,&diffNTRI[0],2,&face_coords[dd],3);
440 A(dd,1) = cblas_ddot(number_nodes,&diffNTRI[1],2,&face_coords[dd],3);
441 }
442 MatrixDouble ATA;
443 ATA = prod(trans(A),A);
444 VectorDouble ATres;
445 ATres = prod(trans(A),res);
446 double det = (ATA(0,0)*ATA(1,1)-ATA(0,1)*ATA(1,0));
447 d_local_coords[0] = -(ATA(1,1)*ATres[0]-ATA(0,1)*ATres[1])/det;
448 d_local_coords[1] = -(ATA(0,0)*ATres[1]-ATA(1,0)*ATres[0])/det;
449 local_coords += d_local_coords;
450 nrm2 = cblas_dnrm2(2,&d_local_coords[0],1);
451 if((it++)>max_it) {
452 cerr << "Warrning: max_it reached in PeriodicNitscheConstrains:: ... ::doAdditionalJobWhenGuassPtsAreCalulated\n";
453 cerr << "nrm2 " << nrm2 << " nrm2_res " << nrm2_res << endl;
454 cerr << "res " << res << endl;
455 cerr << "other_coords " << other_coords << endl;
456 break;
457 }
458 } while(nrm2>eps || nrm2_res>eps);
459 if(local_coords[0]<-eps) {
460 SETERRQ1(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"local_coords[0]=%6.4e<eps",local_coords[0]);
461 }
462 if(local_coords[1]<-eps) {
463 SETERRQ1(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"local_coords[1]=%6.4e<eps",local_coords[1]);
464 }
465 if(local_coords[0]>1+eps) {
466 SETERRQ1(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"local_coords[0]=%6.4e>1+eps",local_coords[0]);
467 }
468 if(local_coords[1]>1+eps) {
469 SETERRQ1(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"local_coords[1]=%6.4e>1+eps",local_coords[1]);
470 }
471 if(local_coords[0]+local_coords[1]>1+eps) {
472 SETERRQ1(
473 PETSC_COMM_SELF,
475 "local_coords[0]+local_coords[1]>1",
476 local_coords[0]+local_coords[1]
477 );
478 }
479 }
480
481 }
482 periodicCommonData.localCoordsMap[other_face].push_back(local_coords);
483 periodicCommonData.inTetFaceGaussPtsNumber[other_face].push_back(gg);
484 }
485 }
488 periodicCommonData.sTress.resize(gg);
491
492 } catch (const std::exception& ex) {
493 ostringstream ss;
494 ss << "throw in method: " << ex.what() << endl;
495 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
496 }
497
498 try {
499 for(
500 map<EntityHandle,vector<int> >::iterator mit = periodicCommonData.inTetFaceGaussPtsNumber.begin();
502 mit++
503 ) {
504 //cerr << mit->first << " " << mit->second.size() << endl;
505 NumeredEntFiniteElement_multiIndex::index<Composite_Name_And_Ent_mi_tag>::type::iterator it;
507 problemPtr->numeredFiniteElements)
508 .get<Composite_Name_And_Ent_mi_tag>()
509 .find(
510 boost::make_tuple(blockData.faceElemName, mit->first));
512 problemPtr->numeredFiniteElements)
513 .get<Composite_Name_And_Ent_mi_tag>()
514 .end()) {
515 SETERRQ1(
516 PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"No finite element found < %s >",
517 blockData.faceElemName.c_str()
518 );
519 }
520 //cerr << *it << endl;
521 boost::shared_ptr<const NumeredEntFiniteElement> face_ptr = *it;
522 periodicFace.copyBasicMethod(*this);
524 periodicFace.nInTheLoop = 0;
525 periodicFace.numeredEntFiniteElementPtr = face_ptr;
526 periodicFace.dataPtr = face_ptr->sPtr->data_dofs;
527 periodicFace.rowPtr = face_ptr->rows_dofs;
528 periodicFace.colPtr = face_ptr->cols_dofs;
529 ierr = periodicFace(); CHKERRQ(ierr);
530
531 Range tets;
532 ierr = mField.get_adjacencies_equality(mit->first,3,tets); CHKERRQ(ierr);
533 if(tets.size()!=1) {
534 SETERRQ1(
535 PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Expected one tet, but is %d",
536 tets.size()
537 );
538 }
540 problemPtr->numeredFiniteElements)
541 .get<Composite_Name_And_Ent_mi_tag>()
542 .find(boost::make_tuple(periodicCommonData.volumeElemName,
543 tets[0]));
545 problemPtr->numeredFiniteElements)
546 .get<Composite_Name_And_Ent_mi_tag>()
547 .end()) {
548 SETERRQ1(
549 PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"No finite element found < %s >",
551 );
552 }
553 boost::shared_ptr<const NumeredEntFiniteElement> tet_ptr = *it;
554 periodicVolume.copyBasicMethod(*this);
556 periodicVolume.nInTheLoop = 0;
557 periodicVolume.numeredEntFiniteElementPtr = tet_ptr;
558 periodicVolume.dataPtr = face_ptr->sPtr->data_dofs;
559 periodicVolume.rowPtr = face_ptr->rows_dofs;
560 periodicVolume.colPtr = face_ptr->cols_dofs;
561 periodicVolume.forFace = mit->first;
562 ierr = periodicVolume(); CHKERRQ(ierr);
563
564 }
565 } catch (const std::exception& ex) {
566 ostringstream ss;
567 ss << "throw in method: " << ex.what() << endl;
568 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
569 }
570
571 PetscFunctionReturn(0);
572 }
573
574 };
575
577
579 vector<Vec> &F;
580
582 string field_name,
583 NitscheMethod::BlockData &nitsche_block_data,
584 NitscheMethod::CommonData &nitsche_common_data,
587 CommonData &periodic_common_data,
588 vector<Vec> &f,
589 bool field_disp = true
590 ):
591 OpCommon(
593 nitsche_block_data,
594 nitsche_common_data,
595 data,
596 common_data,
597 field_disp,
599 ),
600 periodicCommonData(periodic_common_data),
601 F(f) {
602 }
603
604 VectorDouble nF0[6],nF1[6];
605
606
607 VectorDouble voightStrain;
608 VectorDouble dispOnThisSide;
609 VectorDouble dispOnOtherSide;
610
611 MatrixDouble matD0,matD1;
612 static PetscErrorCode calcualteDMatrix(VectorAdaptor coords,MatrixDouble &mat_d) {
613 PetscFunctionBegin;
614 mat_d.resize(3,6,false);
615 mat_d.clear();
616 mat_d(0,0) = coords[0]; mat_d(0,3) = coords[1]*0.5; mat_d(0,5) = coords[2]*0.5;
617 mat_d(1,1) = coords[1]; mat_d(1,3) = coords[0]*0.5; mat_d(1,4) = coords[2]*0.5;
618 mat_d(2,2) = coords[2]; mat_d(2,4) = coords[1]*0.5; mat_d(2,5) = coords[0]*0.5;
619 PetscFunctionReturn(0);
620 }
621
622 PetscErrorCode doWork(
623 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
624 ) {
625 PetscFunctionBegin;
626
627 if(dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) == dAta.tEts.end()) {
628 PetscFunctionReturn(0);
629 }
630
631 if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
632 int nb_dofs = row_data.getIndices().size();
633 const double gamma = nitscheBlockData.gamma;
634 const double phi = nitscheBlockData.phi;
635 double eps = periodicCommonData.eps;
636
637 try {
638
639 int gg = 0;
640 for(int ff = 0;ff<4;ff++) {
641 if(nitscheCommonData.facesFePtr[ff]==NULL) continue;
642 int nb_face_gauss_pts = nitscheCommonData.faceGaussPts[ff].size2();
643 ierr = getFaceRadius(ff); CHKERRQ(ierr);
644 for(int ii = 0;ii<6;ii++) {
645 nF0[ii].resize(nb_dofs,false);
646 nF0[ii].clear();
647 nF1[ii].resize(nb_dofs,false);
648 nF1[ii].clear();
649 }
650 for(int fgg = 0;fgg<nb_face_gauss_pts;fgg++,gg++) {
651 ierr = getGammaH(gamma,gg); CHKERRQ(ierr);
652 double val = getGaussPts()(3,gg);
653 ierr = getJac(row_data,gg,jAc_row); CHKERRQ(ierr);
654 ierr = getTractionVariance(gg,fgg,ff,jAc_row,tRac_v); CHKERRQ(ierr);
655 VectorAdaptor normal = VectorAdaptor(
656 3,ublas::shallow_array_adaptor<double>(
657 3,&nitscheCommonData.faceNormals[ff](fgg,0)
658 )
659 );
660 double area = cblas_dnrm2(3,&normal[0],1);
661
662 voightStrain.resize(6,false);
663 dispOnThisSide.resize(3,false);
664 dispOnOtherSide.resize(3,false);
665 VectorAdaptor x = VectorAdaptor(
666 3,ublas::shallow_array_adaptor<double>(3,&getCoordsAtGaussPts()(gg,0))
667 );
669 VectorDouble other_x = VectorAdaptor(
670 3,ublas::shallow_array_adaptor<double>(3,&periodicCommonData.hoCoordsAtGaussPts[gg][0])
671 );
672 ierr = calcualteDMatrix(other_x,matD1);
673 const double err = 1e-10;
674 double delta = inner_prod(other_x,normal)+inner_prod(x,normal);
675 if(fabs(delta)>err) {
676 cerr << x << endl;
677 cerr << other_x << endl;
678 SETERRQ1(
679 PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"err = %6.4e",delta
680 );
681 }
682
683 for(int ii = 0;ii<6;ii++) {
684 voightStrain.clear();
685 voightStrain[ii] = 1;
686 noalias(dispOnThisSide) = prod(matD0,voightStrain);
687 noalias(dispOnOtherSide) = prod(matD1,voightStrain);
688 for(int dd1 = 0;dd1<nb_dofs/3;dd1++) {
689 double n_val = row_data.getN(gg)[dd1];
690 for(int dd2 = 0;dd2<3;dd2++) {
691 nF0[ii][3*dd1+dd2] -= dispOnThisSide[dd2]*n_val*val*area;
692 nF0[ii][3*dd1+dd2] += (1-eps)*dispOnOtherSide[dd2]*n_val*val*area;
693 }
694 }
695 for(int dd = 0;dd<nb_dofs;dd++) {
696 double dot;
697 dot = cblas_ddot(
698 3,&dispOnThisSide[0],1,&tRac_v(0,dd),tRac_v.size2()
699 );
700 dot -= (1-eps)*cblas_ddot(
701 3,&dispOnOtherSide[0],1,&tRac_v(0,dd),tRac_v.size2()
702 );
703 nF1[ii][dd] += val*phi*dot;
704 }
705 }
706 }
707
708 for(int ii = 0;ii<6;ii++) {
709 nF0[ii] *= -(1/gammaH);
710 ierr = VecSetValues(
711 F[ii],
712 nb_dofs,
713 &row_data.getIndices()[0],
714 &nF0[ii][0],
715 ADD_VALUES
716 ); CHKERRQ(ierr);
717 nF1[ii] *= -1;
718 ierr = VecSetValues(
719 F[ii],
720 nb_dofs,
721 &row_data.getIndices()[0],
722 &nF1[ii][0],
723 ADD_VALUES
724 ); CHKERRQ(ierr);
725 }
726
727 }
728
729 } catch (const std::exception& ex) {
730 ostringstream ss;
731 ss << "throw in method: " << ex.what() << endl;
732 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
733 }
734
735 PetscFunctionReturn(0);
736 }
737
738 };
739
741
743
745 const string field_name,
746 NitscheMethod::BlockData &nitsche_block_data,
747 NitscheMethod::CommonData &nitsche_common_data,
750 CommonData &periodic_common_data,
751 bool field_disp = true
752 ):
753 OpCommon(
755 nitsche_block_data,
756 nitsche_common_data,
757 data,
758 common_data,
759 field_disp,
761 ),
762 periodicCommonData(periodic_common_data) {
763 }
764
765 MatrixDouble kMatrix0,kMatrix1;
766 MatrixDouble kMatrixOff;
767
768 PetscErrorCode doWork(
769 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
770 ) {
771 PetscFunctionBegin;
772
773
774 if(dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) == dAta.tEts.end()) {
775 PetscFunctionReturn(0);
776 }
777 if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
778 int nb_dofs_row = row_data.getIndices().size();
779 const double gamma = nitscheBlockData.gamma;
780 const double phi = nitscheBlockData.phi;
781 double eps = periodicCommonData.eps;
782
783 try {
784
785 int gg = 0;
786 for(int ff = 0;ff<4;ff++) {
787 if(nitscheCommonData.facesFePtr[ff]==NULL) continue;
788 int nb_face_gauss_pts = nitscheCommonData.faceGaussPts[ff].size2();
789 vector<vector<MatrixDouble > > &P = nitscheCommonData.P;
790 P.resize(4);
791 P[ff].resize(nb_face_gauss_pts);
792 for(int fgg = 0;fgg<nb_face_gauss_pts;fgg++) {
793 P[ff][fgg].resize(3,3);
794 P[ff][fgg].clear();
795 for(int dd = 0;dd<3;dd++) {
796 P[ff][fgg](dd,dd) = 1;
797 }
798 }
799 ierr = getFaceRadius(ff); CHKERRQ(ierr);
800 for(int fgg = 0;fgg<nb_face_gauss_pts;fgg++,gg++) {
801 ierr = getGammaH(gamma,gg); CHKERRQ(ierr);
802 double val = getGaussPts()(3,gg);
803 ierr = getJac(row_data,gg,jAc_row); CHKERRQ(ierr);
804 ierr = getTractionVariance(gg,fgg,ff,jAc_row,tRac_v); CHKERRQ(ierr);
805 VectorAdaptor normal = VectorAdaptor(
806 3,ublas::shallow_array_adaptor<double>(3,&nitscheCommonData.faceNormals[ff](fgg,0))
807 );
808 double area = cblas_dnrm2(3,&normal[0],1);
809 const double alpha = 0.5; // mortar parameter
810
811 CommonData::Container::nth_index<0>::type::iterator it,hi_it;
812 it = periodicCommonData.facesContainer.get<0>().lower_bound(gg);
813 hi_it = periodicCommonData.facesContainer.get<0>().upper_bound(gg);
814 for(;it!=hi_it;it++) {
815
816 const ublas::vector<int> &indices = it->iNdices;
817 const VectorDouble &shape = it->shapeFunctions;
818
819 int nb_dofs_col = indices.size();
820 if(indices.size()==0) {
821 continue;
822 }
823 kMatrix0.resize(nb_dofs_row,nb_dofs_col,false);
824 kMatrix0.clear();
825 for(int dd1 = 0;dd1<nb_dofs_row/3;dd1++) {
826 double n_row = row_data.getN()(gg,dd1);
827 for(int dd2 = 0;dd2<nb_dofs_col/3;dd2++) {
828 double n_col = shape[dd2];
829 for(int dd3 = 0;dd3<3;dd3++) {
830 kMatrix0(3*dd1+dd3,3*dd2+dd3) -= (1-eps)*val*n_row*n_col*area;
831 }
832 }
833 }
834 kMatrix1.resize(nb_dofs_row,nb_dofs_col,false);
835 kMatrix1.clear();
836 for(int dd1 = 0;dd1<nb_dofs_row;dd1++) {
837 for(int dd2 = 0;dd2<nb_dofs_col/3;dd2++) {
838 double n_col = shape[dd2];
839 for(int dd3 = 0;dd3<3;dd3++) {
840 kMatrix1(dd1,3*dd2+dd3) += (1-eps)*phi*val*tRac_v(dd3,dd1)*n_col;
841 }
842 }
843 }
844
845 kMatrix0 /= gammaH;
846 ierr = MatSetValues(
847 getFEMethod()->snes_B,
848 nb_dofs_row,
849 &row_data.getIndices()[0],
850 nb_dofs_col,
851 &indices[0],
852 &kMatrix0(0,0),
853 ADD_VALUES
854 ); CHKERRQ(ierr);
855
856 ierr = MatSetValues(
857 getFEMethod()->snes_B,
858 nb_dofs_row,
859 &row_data.getIndices()[0],
860 nb_dofs_col,
861 &indices[0],
862 &kMatrix1(0,0),
863 ADD_VALUES
864 ); CHKERRQ(ierr);
865
866 kMatrix1.resize(nb_dofs_col,nb_dofs_row);
867 kMatrix1.clear();
868 for(int dd1 = 0;dd1<nb_dofs_col/3;dd1++) {
869 double n_col = shape[dd1];
870 for(int dd2 = 0;dd2<3;dd2++) {
871 for(int dd3 = 0;dd3<nb_dofs_row;dd3++) {
872 kMatrix1(3*dd1+dd2,dd3) += (1-eps)*val*n_col*tRac_v(dd2,dd3);
873 }
874 }
875 }
876
877 kMatrix1 *= alpha;
878 ierr = MatSetValues(
879 getFEMethod()->snes_B,
880 nb_dofs_col,
881 &indices[0],
882 nb_dofs_row,
883 &row_data.getIndices()[0],
884 &kMatrix1(0,0),
885 ADD_VALUES
886 ); CHKERRQ(ierr);
887
888 }
889
890 MatrixDouble jac;
891 MatrixDouble trac;
892
893 it = periodicCommonData.volumesContainer.get<0>().lower_bound(gg);
894 hi_it = periodicCommonData.volumesContainer.get<0>().upper_bound(gg);
895 for(;it!=hi_it;it++) {
896
897 int nb = it->iNdices.size();
898 jac.resize(9,nb,false);
899 jac.clear();
900 const MatrixDouble diff_n = it->diffShapeFunctions;
901 const MatrixDouble jac_stress = periodicCommonData.stressJacobian[gg];
902 for(int dd = 0;dd<nb/3;dd++) {
903 for(int rr = 0;rr<3;rr++) {
904 for(int ii = 0;ii<9;ii++) {
905 for(int jj = 0;jj<3;jj++) {
906 jac(ii,3*dd+rr) += jac_stress(ii,3*rr+jj)*diff_n(dd,jj);
907 }
908 }
909 }
910 }
911 trac.resize(3,jac.size2());
912 trac.clear();
913 for(unsigned int dd2 = 0;dd2<jac.size2();dd2++) {
914 for(int nn = 0;nn<3;nn++) {
915 trac(nn,dd2) = -cblas_ddot(3,&jac(3*nn,dd2),jac.size2(),&normal[0],1);
916 }
917 }
918
919 kMatrixOff.resize(nb_dofs_row,nb);
920 kMatrixOff.clear();
921 for(int dd1 = 0;dd1<nb_dofs_row/3;dd1++) {
922 double n_row = row_data.getN()(gg,dd1);
923 for(int dd2 = 0;dd2<3;dd2++) {
924 for(int dd3 = 0;dd3<nb;dd3++) {
925 kMatrixOff(3*dd1+dd2,dd3) += (1-eps)*val*n_row*trac(dd2,dd3);
926 }
927 }
928 }
929
930 kMatrixOff *= (1-alpha);
931 ierr = MatSetValues(
932 getFEMethod()->snes_B,
933 nb_dofs_row,
934 &row_data.getIndices()[0],
935 nb,
936 &it->iNdices[0],
937 &kMatrixOff(0,0),
938 ADD_VALUES
939 ); CHKERRQ(ierr);
940
941 }
942
943 }
944 }
945
946 } catch (const std::exception& ex) {
947 ostringstream ss;
948 ss << "throw in method: " << ex.what() << endl;
949 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
950 }
951
952 PetscFunctionReturn(0);
953 }
954
955 };
956
958
960
962 const string field_name,
965 Vec stress_homo_ghost_vector
966 ):
968 field_name,data,common_data
969 ),
970 stressHomoGhostVector(stress_homo_ghost_vector) {
971 }
972
973 VectorDouble homoStress;
974
975 PetscErrorCode doWork(
976 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
977 ) {
978 PetscFunctionBegin;
979
980
981
982 if(row_type != MBVERTEX) PetscFunctionReturn(0);
983 if(dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) == dAta.tEts.end()) {
984 PetscFunctionReturn(0);
985 }
986
987 try {
988
989 homoStress.resize(6,false);
990 homoStress.clear();
991
992 for(unsigned int gg = 0;gg<row_data.getN().size1();gg++) {
993 const MatrixDouble& stress = commonData.sTress[gg];
994 double val = getGaussPts()(3,gg)*getVolume();
995 if(getHoGaussPtsDetJac().size()>0) {
996 val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
997 }
998 homoStress[0] += stress(0,0)*val;
999 homoStress[1] += stress(1,1)*val;
1000 homoStress[2] += stress(2,2)*val;
1001 homoStress[3] += (stress(0,1)+stress(1,0))*0.5*val;
1002 homoStress[4] += (stress(1,2)+stress(2,1))*0.5*val;
1003 homoStress[5] += (stress(0,2)+stress(2,0))*0.5*val;
1004 }
1005
1006 const int indices_6[6] = {0, 1, 2, 3, 4, 5};
1007 ierr = VecSetValues(
1008 stressHomoGhostVector,6,indices_6,&homoStress[0],ADD_VALUES
1009 ); CHKERRQ(ierr);
1010
1011 } catch (const std::exception& ex) {
1012 ostringstream ss;
1013 ss << "throw in method: " << ex.what() << endl;
1014 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1015 }
1016
1017 PetscFunctionReturn(0);
1018 }
1019
1020 };
1021
1023
1025
1027 const string field_name,
1028 NitscheMethod::BlockData &nitsche_block_data,
1029 NitscheMethod::CommonData &nitsche_common_data,
1032 Vec stress_homo_ghost_vector,
1033 bool field_disp = 0
1034 ):
1036 field_name,
1037 nitsche_block_data,
1038 nitsche_common_data,
1039 data,
1040 common_data,
1041 field_disp,
1043 ),
1044 stressHomoGhostVector(stress_homo_ghost_vector) {
1045 }
1046
1047 MatrixDouble matD;
1048 VectorDouble tRaction;
1049 VectorDouble homoStress;
1050
1051 VectorDouble resDisp;
1052 VectorDouble dispOnOtherSide;
1053
1054 PetscErrorCode doWork(
1055 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
1056 ) {
1057 PetscFunctionBegin;
1058
1059
1060
1061 if(row_type != MBVERTEX) PetscFunctionReturn(0);
1062 if(dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) == dAta.tEts.end()) {
1063 PetscFunctionReturn(0);
1064 }
1065
1066 try {
1067
1068 homoStress.resize(6,false);
1069 homoStress.clear();
1070
1071 int gg = 0;
1072 for(int ff = 0;ff<4;ff++) {
1073 if(nitscheCommonData.facesFePtr[ff]==NULL) continue;
1074 int nb_face_gauss_pts = nitscheCommonData.faceGaussPts[ff].size2();
1075 for(int fgg = 0;fgg<nb_face_gauss_pts;fgg++,gg++) {
1076 double val = getGaussPts()(3,gg);
1077 const MatrixDouble& stress = commonData.sTress[gg];
1078 VectorAdaptor normal = VectorAdaptor(
1079 3,ublas::shallow_array_adaptor<double>(
1080 3,&nitscheCommonData.faceNormals[ff](fgg,0)
1081 )
1082 );
1083
1084 VectorAdaptor x = VectorAdaptor(
1085 3,ublas::shallow_array_adaptor<double>(
1086 3,&getCoordsAtGaussPts()(gg,0)
1087 )
1088 );
1089
1090 tRaction.resize(3,false);
1091 noalias(tRaction) = prod(normal,stress);
1092
1093 homoStress[0] += x[0]*tRaction[0]*val;
1094 homoStress[1] += x[1]*tRaction[1]*val;
1095 homoStress[2] += x[2]*tRaction[2]*val;
1096 homoStress[3] += (x[0]*tRaction[1]+x[1]*tRaction[0])*0.5*val;
1097 homoStress[4] += (x[1]*tRaction[2]+x[2]*tRaction[1])*0.5*val;
1098 homoStress[5] += (x[0]*tRaction[2]+x[2]*tRaction[0])*0.5*val;
1099 }
1100 }
1101
1102 const int indices_6[6] = {0, 1, 2, 3, 4, 5};
1103 ierr = VecSetValues(
1104 stressHomoGhostVector,6,indices_6,&homoStress[0],ADD_VALUES
1105 ); CHKERRQ(ierr);
1106
1107 if(gg != (int)getGaussPts().size2()) {
1108 SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong number of gauss pts");
1109 }
1110
1111
1112 } catch (const std::exception& ex) {
1113 ostringstream ss;
1114 ss << "throw in method: " << ex.what() << endl;
1115 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1116 }
1117
1118 PetscFunctionReturn(0);
1119 }
1120
1121 };
1122
1123
1124};
1125
1126#endif //__NITSCHE_PERIODIC_METHOD_HPP__
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
NumeredEntFiniteElement_multiIndex & returnNumeredEntFiniteElement_multiIndex(T &fes)
static PetscErrorCode ierr
static const double eps
#define CHKERRQ_MOAB(a)
check error code of MoAB function
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
constexpr int order
PetscErrorCode ShapeDiffMBTRIQ(double *diffN, const double *X, const double *Y, const int G_DIM)
Definition fem_tools.c:795
PetscErrorCode ShapeFaceNormalMBTRI(double *diffN, const double *coords, double *normal)
Definition fem_tools.c:229
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition fem_tools.c:194
PetscErrorCode ShapeMBTRIQ(double *N, const double *X, const double *Y, const int G_DIM)
Definition fem_tools.c:780
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition fem_tools.c:182
static double phi
static const double face_coords[4][9]
double tol
constexpr AssemblyType A
constexpr auto field_name
static constexpr double delta
const int N
Definition speed_test.cpp:3
const Problem * problemPtr
raw pointer to problem
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
@ OPROW
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
Block data for Nitsche method.
double gamma
Penalty term, see .
string faceElemName
name of element face
double phi
Nitsche method parameter, see .
Common data shared between finite element operators.
std::vector< boost::shared_ptr< const NumeredEntFiniteElement > > facesFePtr
std::vector< MatrixDouble > faceNormals
std::vector< VectorDouble > rAy
std::vector< std::vector< MatrixDouble > > P
projection matrix
std::vector< MatrixDouble > hoCoordsAtGaussPts
std::vector< MatrixDouble > faceGaussPts
Definition of volume element.
MyVolumeFE(MoFEM::Interface &m_field, BlockData &block_data, CommonData &common_data)
virtual MoFEMErrorCode getFaceRadius(int ff)
virtual MoFEMErrorCode getGammaH(double gamma, int gg)
Calculate jacobian and variation of tractions.
NonlinearElasticElement::CommonData & commonData
OpCommon(const std::string field_name, BlockData &nitsche_block_data, CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, bool field_disp, const char type)
NonlinearElasticElement::BlockData & dAta
MoFEMErrorCode getTractionVariance(int gg, int fgg, int ff, MatrixDouble &jac, MatrixDouble &trac)
Basic implementation of Nitsche's method.
data for calculation heat conductivity and heat capacity elements
Range tEts
constrains elements in block set
common data used by volume elements
std::vector< MatrixDouble > jacStress
this is simply material tangent operator
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
structure grouping operators and data used for calculation of nonlinear elastic element
map< EntityHandle, vector< VectorDouble > > localCoordsMap
map< EntityHandle, vector< int > > inTetFaceGaussPtsNumber
multi_index_container< MultiIndexData, indexed_by< ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData, int, MultiIndexData::gG) >, ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData, int, MultiIndexData::sIde) >, ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData, EntityType, MultiIndexData::tYpe) >, ordered_unique< composite_key< MultiIndexData, member< MultiIndexData, int,&MultiIndexData::gG >, member< MultiIndexData, int,&MultiIndexData::sIde >, member< MultiIndexData, EntityType,&MultiIndexData::tYpe > > > > > Container
map< EntityHandle, vector< int > > inTetTetGaussPtsNumber
MyNitscheVolume(MoFEM::Interface &m_field, NitscheMethod::BlockData &block_data, NitscheMethod::CommonData &common_data, CommonData &periodic_common_data)
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpCalculateHomogenisedStressSurfaceIntegral(const string field_name, NitscheMethod::BlockData &nitsche_block_data, NitscheMethod::CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, Vec stress_homo_ghost_vector, bool field_disp=0)
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpCalculateHomogenisedStressVolumeIntegral(const string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, Vec stress_homo_ghost_vector)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetFaceData(CommonData &common_data, bool field_disp=true)
NonlinearElasticElement::CommonData & commonData
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetVolumeData(NonlinearElasticElement::CommonData &common_data, CommonData &periodic_common_data)
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpLhsPeriodicNormal(const string field_name, NitscheMethod::BlockData &nitsche_block_data, NitscheMethod::CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, CommonData &periodic_common_data, bool field_disp=true)
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpRhsPeriodicNormal(string field_name, NitscheMethod::BlockData &nitsche_block_data, NitscheMethod::CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, CommonData &periodic_common_data, vector< Vec > &f, bool field_disp=true)
static PetscErrorCode calcualteDMatrix(VectorAdaptor coords, MatrixDouble &mat_d)
PeriodicFace(MoFEM::Interface &m_field, CommonData &common_data)
PeriodicVolume(MoFEM::Interface &m_field, CommonData &common_data)