v0.14.0
Loading...
Searching...
No Matches
NitscheMethod.hpp
Go to the documentation of this file.
1/** \file NitscheMethod.hpp
2 * \ingroup nitsche_method
3 * \brief Basic implementation of Nitsche method
4 *
5 * \note This is quite an old implementation, kept here for compatibility with
6 * older users modules. Currently, MoFEM has a generic implementation for
7 * integration over the skeleton, and that one should be used, instead one
8 * presented below.
9 *
10 */
11
12/*
13 * This file is part of MoFEM.
14 * MoFEM is free software: you can redistribute it and/or modify it under
15 * the terms of the GNU Lesser General Public License as published by the
16 * Free Software Foundation, either version 3 of the License, or (at your
17 * option) any later version.
18 *
19 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
20 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
22 * License for more details.
23 *
24 * You should have received a copy of the GNU Lesser General Public
25 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
26
27#ifndef __NITCHE_BOUNDARY_CONDITIONS_HPP__
28#define __NITCHE_BOUNDARY_CONDITIONS_HPP__
29
30// This is fix to have back compatibility
31template <typename T>
34 return fes;
35}
36
37template <>
41 return fes;
42}
43
44template <>
47 boost::shared_ptr<NumeredEntFiniteElement_multiIndex>>(
48 boost::shared_ptr<NumeredEntFiniteElement_multiIndex> &fes) {
49 return *fes;
50}
51
52/** \brief Basic implementation of Nitsche's method
53 * \ingroup nitsche_method
54
55 This implementation have been used to enforce periodic constrains, see \ref
56 nitsche_periodic. For theoretical foundations of the method explanations see
57 \cite nitsche_method_hal and \cite juntunen2009nitsche. For scratch book
58 derivations see <a href="nitsche_bc_for_ch.pdf" target="_blank"><b>link</b></a>.
59
60*/
62
67 int getRule(int order) { return order + addToRule; }
68 /*int getRule(int order) { return -1; }
69 MoFEMErrorCode setGaussPts(int order) {
70 MoFEMFunctionBeginHot;
71 int rule = order+addToRule+1;
72 int nb_gauss_pts = triangle_ncc_order_num(rule);
73 gaussPts.resize(3,nb_gauss_pts,false);
74 triangle_ncc_rule(rule,nb_gauss_pts,&gaussPts(0,0),&gaussPts(2,0));
75 MoFEMFunctionReturnHot(0);
76 }*/
77 };
78
79 /** \brief Block data for Nitsche method
80 * \ingroup nitsche_method
81 */
82 struct BlockData {
83 double gamma; ///< Penalty term, see \cite nitsche_method_hal
84 double phi; ///< Nitsche method parameter, see \cite nitsche_method_hal
85 string faceElemName; ///< name of element face
86 Range fAces; ///< faces on which constrain is applied
87 };
88
89 /** \brief Common data shared between finite element operators
90 */
91 struct CommonData {
93 std::vector<EntityHandle> fAces;
94 std::vector<boost::shared_ptr<const NumeredEntFiniteElement>> facesFePtr;
95 std::vector<VectorDouble> cOords;
96 std::vector<MatrixDouble> faceNormals;
97 std::vector<MatrixDouble> faceGaussPts;
98 std::vector<std::vector<int>> inTetFaceGaussPtsNumber;
99 std::vector<MatrixDouble> coordsAtGaussPts;
100 std::vector<MatrixDouble> hoCoordsAtGaussPts;
101 std::vector<VectorDouble> rAy;
103
104 /** \brief projection matrix
105
106 First index is face index, then integration point and
107 projection matrix 3x3
108
109 */
110 std::vector<std::vector<MatrixDouble>> P;
111
112 /** \brief derivative of projection matrix in respect DoFs
113 This is EntityType, face, gauss point at face.
114 Matrix has 3 rows (components of displacements)
115 Matrix has columns equal to number of DoFs on entity
116 */
117 std::map<EntityType, std::vector<std::vector<MatrixDouble>>> dP;
118
120
122 int gG;
123 int sIde;
125 ublas::vector<int> iNdices;
126 ublas::vector<int> dofOrders;
127 VectorDouble shapeFunctions;
128 MatrixDouble diffShapeFunctions;
129 MultiIndexData(int gg, int side, EntityType type)
130 : gG(gg), sIde(side), tYpe(type) {}
131 };
132
133 typedef multi_index_container<
134 MultiIndexData,
135 indexed_by<
136 ordered_non_unique<BOOST_MULTI_INDEX_MEMBER(MultiIndexData, int,
138 ordered_non_unique<BOOST_MULTI_INDEX_MEMBER(MultiIndexData, int,
140 ordered_non_unique<BOOST_MULTI_INDEX_MEMBER(
141 MultiIndexData, EntityType, MultiIndexData::tYpe)>,
142 ordered_unique<composite_key<
143 MultiIndexData,
144 member<MultiIndexData, int, &MultiIndexData::gG>,
145 member<MultiIndexData, int, &MultiIndexData::sIde>,
146 member<MultiIndexData, EntityType, &MultiIndexData::tYpe>>>>>
149 };
150
151 /** \brief Get integration pts data on face
152 */
155
157
160 "DISPLACEMENT", OPROW),
161 commonData(common_data) {}
162
163 MoFEMErrorCode doWork(int side, EntityType type,
164 DataForcesAndSourcesCore::EntData &data) {
166
167 int faceInRespectToTet = getFEMethod()->nInTheLoop;
168 int nb_face_gauss_pts = getGaussPts().size2();
169
170 try {
171 if (type == MBVERTEX) {
172 commonData.faceGaussPts.resize(4);
173 commonData.faceGaussPts[faceInRespectToTet] = getGaussPts();
174 for (int fgg = 0; fgg < nb_face_gauss_pts; fgg++) {
175 int gg = commonData.nbTetGaussPts;
176 commonData.inTetFaceGaussPtsNumber[faceInRespectToTet].push_back(
177 gg);
179 }
180 commonData.faceNormals.resize(4);
181 commonData.faceNormals[faceInRespectToTet] =
182 0.5 * getNormalsAtGaussPtss();
184 commonData.hoCoordsAtGaussPts[faceInRespectToTet] =
186 commonData.cOords.resize(4);
187 commonData.cOords[faceInRespectToTet] = getCoords();
189 commonData.coordsAtGaussPts[faceInRespectToTet] =
191 commonData.rAy.resize(4);
192 commonData.rAy[faceInRespectToTet] = -getNormal();
193 commonData.rAy[faceInRespectToTet] /= norm_2(getNormal());
194 }
195 } catch (const std::exception &ex) {
196 std::ostringstream ss;
197 ss << "throw in method: " << ex.what() << std::endl;
198 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
199 }
200
201 try {
202 for (int fgg = 0; fgg < nb_face_gauss_pts; fgg++) {
203 int gg = commonData.inTetFaceGaussPtsNumber[faceInRespectToTet][fgg];
204 // std::cerr << fgg << " " << gg << " " << side << " " << type <<
205 // std::endl;
206 CommonData::MultiIndexData gauss_pt_data(gg, side, type);
207 std::pair<CommonData::Container::iterator, bool> p;
208 p = commonData.facesContainer.insert(gauss_pt_data);
209 if (!p.second) {
210 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
211 "data not inserted");
212 }
214 const_cast<CommonData::MultiIndexData &>(*p.first);
215 VectorDouble &shape_fun = p_data.shapeFunctions;
216 int nb_shape_fun = data.getN().size2();
217 shape_fun.resize(nb_shape_fun);
218 // std::cerr << "nb_shape_fun " << nb_shape_fun << std::endl;
219 if (nb_shape_fun) {
220 cblas_dcopy(nb_shape_fun, &data.getN()(fgg, 0), 1, &shape_fun[0],
221 1);
222 }
223 p_data.iNdices = data.getIndices();
224 p_data.dofOrders.resize(data.getFieldDofs().size(), false);
225 for (unsigned int dd = 0; dd < data.getFieldDofs().size(); dd++) {
226 p_data.dofOrders[dd] = data.getFieldDofs()[dd]->getDofOrder();
227 }
228 // std::cerr << shape_fun << std::endl;
229 }
230 } catch (const std::exception &ex) {
231 std::ostringstream ss;
232 ss << "throw in method: " << ex.what() << std::endl;
233 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
234 }
236 }
237 };
238
239 /// \brief Definition of volume element
241
246
250 }
251
252 MyVolumeFE(MoFEM::Interface &m_field, BlockData &block_data,
253 CommonData &common_data)
255 blockData(block_data), commonData(common_data), faceFE(m_field),
256 addToRule(0) {
258 }
259
260 int getRule(int order) { return -1; };
261
262 MoFEMErrorCode setGaussPts(int order) {
263
265
266 gaussPts.resize(4, 0, false);
267
268 try {
270 commonData.fAces.resize(4);
272 for (int ff = 0; ff < 4; ff++) {
273 EntityHandle face;
274 rval = mField.get_moab().side_element(tet, 2, ff, face);
275 CHKERRG(rval);
276 if (blockData.fAces.find(face) != blockData.fAces.end()) {
277 commonData.fAces[ff] = face;
279 } else {
280 commonData.fAces[ff] = 0;
281 }
282 }
283 } catch (const std::exception &ex) {
284 std::ostringstream ss;
285 ss << "throw in method: " << ex.what() << std::endl;
286 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
287 }
288
289 try {
290 commonData.facesFePtr.resize(4);
291 for (int ff = 0; ff < 4; ff++) {
292 if (commonData.fAces[ff] != 0) {
293 NumeredEntFiniteElement_multiIndex::index<
294 Composite_Name_And_Ent_mi_tag>::type::iterator it,
295 hi_it;
296
298 problemPtr->numeredFiniteElements)
299 .get<Composite_Name_And_Ent_mi_tag>()
300 .lower_bound(boost::make_tuple(blockData.faceElemName,
301 commonData.fAces[ff]));
303 problemPtr->numeredFiniteElements)
304 .get<Composite_Name_And_Ent_mi_tag>()
305 .upper_bound(boost::make_tuple(blockData.faceElemName,
306 commonData.fAces[ff]));
308 problemPtr->numeredFiniteElements)
309 .get<Composite_Name_And_Ent_mi_tag>()
310 .end()) {
311 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
312 "No finite element found < %s >",
313 blockData.faceElemName.c_str());
314 }
315 commonData.facesFePtr[ff] = *it;
316 } else {
317 commonData.facesFePtr[ff].reset();
318 }
319 }
324 for (int ff = 0; ff < 4; ff++) {
325 if (commonData.facesFePtr[ff]) {
326 boost::shared_ptr<const NumeredEntFiniteElement> faceFEPtr =
328 faceFE.copyBasicMethod(*this);
330 faceFE.nInTheLoop = ff;
332 faceFE.dataPtr = faceFEPtr->sPtr->data_dofs;
333 faceFE.rowPtr = faceFEPtr->rows_dofs;
334 faceFE.colPtr = faceFEPtr->cols_dofs;
336 ierr = faceFE();
337 CHKERRG(ierr);
338 }
339 }
340 } catch (const std::exception &ex) {
341 std::ostringstream ss;
342 ss << "throw in method: " << ex.what() << std::endl;
343 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
344 }
345
346 try {
347 int nb_gauss_pts = 0;
348 for (int ff = 0; ff < 4; ff++) {
349 if (!commonData.facesFePtr[ff])
350 continue;
351 nb_gauss_pts += commonData.faceGaussPts[ff].size2();
352 }
353 gaussPts.resize(4, nb_gauss_pts, false);
354 const double coords_tet[12] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
355 int gg = 0;
356 for (int ff = 0; ff < 4; ff++) {
357 if (!commonData.facesFePtr[ff])
358 continue;
359 int nb_gauss_face_pts = commonData.faceGaussPts[ff].size2();
360 // std::cerr << "nb_gauss_face_pts " << nb_gauss_face_pts <<
361 // std::endl;
362 for (int fgg = 0; fgg < nb_gauss_face_pts; fgg++, gg++) {
363 // std::cerr << ff << " gg " << gg << " fgg " << fgg << std::endl;
364 CommonData::Container::nth_index<3>::type::iterator sit;
365 sit = commonData.facesContainer.get<3>().find(
366 boost::make_tuple(gg, 0, MBVERTEX));
367 const VectorDouble &shape_fun = sit->shapeFunctions;
368 // std::cerr << shape_fun << std::endl;
369
370 for (int dd = 0; dd < 3; dd++) {
371 gaussPts(dd, gg) =
372 shape_fun[0] * coords_tet[3 * dataH1.facesNodes(ff, 0) + dd] +
373 shape_fun[1] * coords_tet[3 * dataH1.facesNodes(ff, 1) + dd] +
374 shape_fun[2] * coords_tet[3 * dataH1.facesNodes(ff, 2) + dd];
375 }
376 gaussPts(3, gg) = commonData.faceGaussPts[ff](2, fgg);
377 }
378 }
379 } catch (const std::exception &ex) {
380 std::ostringstream ss;
381 ss << "throw in method: " << ex.what() << std::endl;
382 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
383 }
384
386 CHKERRG(ierr);
387
389 }
390 };
391
392 /** \brief Basic operated shared between all Nitsche operators
393 */
396
400
401 OpBasicCommon(const std::string field_name, BlockData &nitsche_block_data,
402 CommonData &nitsche_common_data, bool field_disp,
403 const char type)
404 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,
405 type),
406 nitscheBlockData(nitsche_block_data),
407 nitscheCommonData(nitsche_common_data), fieldDisp(field_disp) {
408 sYmm = false;
409 }
410
412 virtual MoFEMErrorCode getFaceRadius(int ff) {
414 VectorDouble &coords = nitscheCommonData.cOords[ff];
415 double center[3];
416 tricircumcenter3d_tp(&coords[0], &coords[3], &coords[6], center, NULL,
417 NULL);
418 cblas_daxpy(3, -1, &coords[0], 1, center, 1);
419 faceRadius = cblas_dnrm2(3, center, 1);
421 }
422 double gammaH;
423 virtual MoFEMErrorCode getGammaH(double gamma, int gg) {
425 gammaH = gamma;
426 // gammaH*= faceRadius;
428 }
429 };
430
431 /** \brief Calculate jacobian and variation of tractions
432 */
433 struct OpCommon : public OpBasicCommon {
434
437
438 OpCommon(const std::string field_name, BlockData &nitsche_block_data,
439 CommonData &nitsche_common_data,
441 NonlinearElasticElement::CommonData &common_data, bool field_disp,
442 const char type)
443 : OpBasicCommon(field_name, nitsche_block_data, nitsche_common_data,
444 field_disp, type),
445 dAta(data), commonData(common_data) {}
446
447 VectorDouble dIsp;
448 VectorDouble tRaction;
449 MatrixDouble jAc_row;
450 MatrixDouble jAc_col;
451 MatrixDouble tRac_v;
452 MatrixDouble tRac_u;
453
454 virtual MoFEMErrorCode calculateP(int gg, int fgg, int ff) {
457 }
458
459 MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &data, int gg,
460 MatrixDouble &jac) {
462 try {
463 int nb = data.getFieldData().size();
464 jac.resize(9, nb, false);
465 jac.clear();
466 const MatrixAdaptor diffN = data.getDiffN(gg, nb / 3);
467 MatrixDouble &jac_stress = commonData.jacStress[gg];
468 for (int dd = 0; dd < nb / 3; dd++) {
469 for (int rr = 0; rr < 3; rr++) {
470 for (int ii = 0; ii < 9; ii++) {
471 for (int jj = 0; jj < 3; jj++) {
472 jac(ii, 3 * dd + rr) +=
473 jac_stress(ii, 3 * rr + jj) * diffN(dd, jj);
474 }
475 }
476 }
477 }
478 } catch (const std::exception &ex) {
479 std::ostringstream ss;
480 ss << "throw in method: " << ex.what() << std::endl;
481 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
482 }
484 }
485
486 MoFEMErrorCode getTractionVariance(int gg, int fgg, int ff,
487 MatrixDouble &jac, MatrixDouble &trac) {
489 try {
490 VectorAdaptor normal = VectorAdaptor(
491 3, ublas::shallow_array_adaptor<double>(
492 3, &nitscheCommonData.faceNormals[ff](fgg, 0)));
493 trac.resize(3, jac.size2());
494 trac.clear();
495 for (unsigned int dd2 = 0; dd2 < jac.size2(); dd2++) {
496 for (int nn = 0; nn < 3; nn++) {
497 trac(nn, dd2) =
498 cblas_ddot(3, &jac(3 * nn, dd2), jac.size2(), &normal[0], 1);
499 }
500 }
501 } catch (const std::exception &ex) {
502 std::ostringstream ss;
503 ss << "throw in method: " << ex.what() << std::endl;
504 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
505 }
507 }
508 };
509
510 /** \brief Calculate Nitsche method terms on left hand side
511 * \ingroup nitsche_method
512 */
513 struct OpLhsNormal : public OpCommon {
514
515 OpLhsNormal(const std::string field_name, BlockData &nitsche_block_data,
516 CommonData &nitsche_common_data,
519 bool field_disp)
520 : OpCommon(field_name, nitsche_block_data, nitsche_common_data, data,
521 common_data, field_disp, UserDataOperator::OPROWCOL) {}
522
523 MatrixDouble kMatrix, kMatrix0, kMatrix1;
524 std::vector<MatrixDouble> kMatrixFace, kMatrixFace0, kMatrixFace1;
525
526 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
527 EntityType col_type,
528 DataForcesAndSourcesCore::EntData &row_data,
529 DataForcesAndSourcesCore::EntData &col_data) {
531
532 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
533 dAta.tEts.end()) {
535 }
536 if (row_data.getIndices().size() == 0)
538 if (col_data.getIndices().size() == 0)
540 int nb_dofs_row = row_data.getIndices().size();
541 int nb_dofs_col = col_data.getIndices().size();
542 double gamma = nitscheBlockData.gamma;
543 double phi = nitscheBlockData.phi;
544
545 kMatrix0.resize(nb_dofs_row, nb_dofs_col, false);
546 kMatrix1.resize(nb_dofs_row, nb_dofs_col, false);
547 kMatrix.resize(nb_dofs_row, nb_dofs_col, false);
548 kMatrix.clear();
549
550 try {
551
552 int gg = 0;
553 for (int ff = 0; ff < 4; ff++) {
555 continue;
556 int nb_face_gauss_pts = nitscheCommonData.faceGaussPts[ff].size2();
557 ierr = getFaceRadius(ff);
558 CHKERRG(ierr);
559 kMatrix0.clear();
560 kMatrix1.clear();
561 for (int fgg = 0; fgg < nb_face_gauss_pts; fgg++, gg++) {
562 ierr = getGammaH(gamma, gg);
563 CHKERRG(ierr);
564 double val = getGaussPts()(3, gg);
565 ierr = getJac(row_data, gg, jAc_row);
566 CHKERRG(ierr);
567 ierr = getTractionVariance(gg, fgg, ff, jAc_row, tRac_v);
568 CHKERRG(ierr);
569 ierr = getJac(col_data, gg, jAc_col);
570 CHKERRG(ierr);
571 ierr = getTractionVariance(gg, fgg, ff, jAc_col, tRac_u);
572 CHKERRG(ierr);
573 VectorAdaptor normal = VectorAdaptor(
574 3, ublas::shallow_array_adaptor<double>(
575 3, &nitscheCommonData.faceNormals[ff](fgg, 0)));
576 double area = cblas_dnrm2(3, &normal[0], 1);
577
578 ierr = calculateP(gg, fgg, ff);
579 CHKERRG(ierr);
580 MatrixDouble &P = nitscheCommonData.P[ff][fgg];
581
582 // P
583 for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
584 double n_row = row_data.getN()(gg, dd1);
585 for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
586 double n_col = col_data.getN()(gg, dd2);
587 for (int dd3 = 0; dd3 < 3; dd3++) {
588 for (int dd4 = 0; dd4 < 3; dd4++) {
589 if (!P(dd3, dd4))
590 continue;
591 kMatrix0(3 * dd1 + dd3, 3 * dd2 + dd4) +=
592 val * n_row * P(dd3, dd4) * n_col * area;
593 }
594 }
595 }
596 }
597 for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
598 double n_row = row_data.getN()(gg, dd1);
599 for (int dd2 = 0; dd2 < nb_dofs_col; dd2++) {
600 for (int dd3 = 0; dd3 < 3; dd3++) {
601 double t = cblas_ddot(3, &P(dd3, 0), 1, &tRac_u(0, dd2),
602 tRac_u.size2());
603 kMatrix1(3 * dd1 + dd3, dd2) -= val * n_row * t;
604 }
605 }
606 }
607 for (int dd1 = 0; dd1 < nb_dofs_row; dd1++) {
608 for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
609 double n_col = col_data.getN()(gg, dd2);
610 for (int dd3 = 0; dd3 < 3; dd3++) {
611 double t = cblas_ddot(3, &P(0, dd3), 3, &tRac_v(0, dd1),
612 tRac_v.size2());
613 kMatrix1(dd1, 3 * dd2 + dd3) -= phi * val * t * n_col;
614 }
615 }
616 }
617 // dP
618 if (nitscheCommonData.dP[col_type].empty() != 0) {
619 if (nitscheCommonData.dP[col_type].size() == 4) {
620 if (nitscheCommonData.dP[col_type][ff].size() ==
621 (unsigned int)nb_face_gauss_pts) {
622 MatrixDouble &dP = nitscheCommonData.dP[col_type][ff][fgg];
623 if (dP.size1() == 3 &&
624 dP.size2() == (unsigned int)nb_dofs_col) {
625 VectorDouble &u =
628 for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
629 double n_row = row_data.getN()(gg, dd1);
630 for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
631 for (int dd3 = 0; dd3 < 3; dd3++) {
632 for (int dd4 = 0; dd4 < 3; dd4++) {
633 kMatrix0(3 * dd1 + dd3, 3 * dd2 + dd4) +=
634 val * n_row * u[dd4] * dP(dd4, 3 * dd2 + dd4);
635 }
636 }
637 }
638 }
639 VectorDouble t =
640 prod(trans(commonData.sTress[gg]),
641 normal); // this is Piola-Kirchhoff I stress
642 for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
643 double n_row = row_data.getN()(gg, dd1);
644 for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
645 for (int dd3 = 0; dd3 < 3; dd3++) {
646 for (int dd4 = 0; dd4 < 3; dd4++) {
647 kMatrix1(3 * dd1 + dd3, 3 * dd2 + dd4) +=
648 val * n_row * t[dd4] * dP(dd4, 3 * dd2 + dd4);
649 }
650 }
651 }
652 }
653 }
654 }
655 }
656 }
657 }
658
659 kMatrix0 /= gammaH;
660 noalias(kMatrix) += kMatrix0 + kMatrix1;
661 }
662
663 if (gg != (int)getGaussPts().size2()) {
664 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
665 "wrong number of gauss pts");
666 }
667
668 ierr = MatSetValues(
669 getFEMethod()->snes_B, nb_dofs_row, &row_data.getIndices()[0],
670 nb_dofs_col, &col_data.getIndices()[0], &kMatrix(0, 0), ADD_VALUES);
671 CHKERRG(ierr);
672
673 } catch (const std::exception &ex) {
674 std::ostringstream ss;
675 ss << "throw in method: " << ex.what() << std::endl;
676 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
677 }
678
680 }
681 };
682
683 /** \brief Calculate Nitsche method terms on right hand side
684 * \ingroup nitsche_method
685 */
686 struct OpRhsNormal : public OpCommon {
687
688 OpRhsNormal(const std::string field_name, BlockData &nitsche_block_data,
689 CommonData &nitsche_common_data,
692 bool field_disp)
693 : OpCommon(field_name, nitsche_block_data, nitsche_common_data, data,
694 common_data, field_disp, UserDataOperator::OPROW) {}
695
696 VectorDouble nF;
697
698 MoFEMErrorCode doWork(int row_side, EntityType row_type,
699 DataForcesAndSourcesCore::EntData &row_data) {
701
702 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
703 dAta.tEts.end()) {
705 }
706 if (row_data.getIndices().size() == 0)
708 int nb_dofs_row = row_data.getIndices().size();
709 double gamma = nitscheBlockData.gamma;
710 double phi = nitscheBlockData.phi;
711
712 try {
713
714 nF.resize(nb_dofs_row, false);
715 nF.clear();
716
717 int gg = 0;
718 for (int ff = 0; ff < 4; ff++) {
720 continue;
721 int nb_face_gauss_pts = nitscheCommonData.faceGaussPts[ff].size2();
722 ierr = getFaceRadius(ff);
723 CHKERRG(ierr);
724 for (int fgg = 0; fgg < nb_face_gauss_pts; fgg++, gg++) {
725
726 ierr = getGammaH(gamma, gg);
727 CHKERRG(ierr);
728 double val = getGaussPts()(3, gg);
729 ierr = getJac(row_data, gg, jAc_row);
730 CHKERRG(ierr);
731 ierr = getTractionVariance(gg, fgg, ff, jAc_row, tRac_v);
732 CHKERRG(ierr);
733 VectorAdaptor normal = VectorAdaptor(
734 3, ublas::shallow_array_adaptor<double>(
735 3, &nitscheCommonData.faceNormals[ff](fgg, 0)));
736 double area = cblas_dnrm2(3, &normal[0], 1);
737
738 ierr = calculateP(gg, fgg, ff);
739 CHKERRG(ierr);
740 MatrixDouble &P = nitscheCommonData.P[ff][fgg];
741
742 VectorDouble &u =
744 VectorDouble t = prod(trans(commonData.sTress[gg]),
745 normal); // this is Piola-Kirchhoff I stress
746
747 for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
748 double n_row = row_data.getN()(gg, dd1);
749 for (int dd2 = 0; dd2 < 3; dd2++) {
750 nF[3 * dd1 + dd2] +=
751 gammaH *
752 (val * area * n_row *
753 (P(dd2, 0) * u[0] + P(dd2, 1) * u[1] + P(dd2, 2) * u[2]));
754 nF[3 * dd1 + dd2] +=
755 val * n_row *
756 (P(dd2, 0) * t[0] + P(dd2, 1) * t[1] + P(dd2, 2) * t[2]);
757 nF[3 * dd1 + dd2] +=
758 phi * tRac_v(dd2, 3 * dd1 + dd2) *
759 (P(dd2, 0) * u[0] + P(dd2, 1) * u[1] + P(dd2, 2) * u[2]);
760 }
761 }
762 }
763 }
764
765 ierr = VecSetValues(getFEMethod()->snes_f, nb_dofs_row,
766 &row_data.getIndices()[0], &nF[0], ADD_VALUES);
767 CHKERRG(ierr);
768
769 } catch (const std::exception &ex) {
770 std::ostringstream ss;
771 ss << "throw in method: " << ex.what() << std::endl;
772 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
773 }
774
776 }
777 };
778};
779
780#endif // __NITCHE_BOUNDARY_CONDITIONS_HPP__
781
782/***************************************************************************/ /**
783* \defgroup nitsche_method Nitsche Method
784* \ingroup user_modules
785******************************************************************************/
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
NumeredEntFiniteElement_multiIndex & returnNumeredEntFiniteElement_multiIndex< NumeredEntFiniteElement_multiIndex >(NumeredEntFiniteElement_multiIndex &fes)
NumeredEntFiniteElement_multiIndex & returnNumeredEntFiniteElement_multiIndex(T &fes)
static PetscErrorCode ierr
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
static double phi
multi_index_container< boost::shared_ptr< NumeredEntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_EntFiniteElement, UId, &NumeredEntFiniteElement::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > > > > > NumeredEntFiniteElement_multiIndex
MultiIndex for entities for NumeredEntFiniteElement.
const double T
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
MoFEMErrorCode copyBasicMethod(const BasicMethod &basic)
Copy data from other base method to this base method.
int nInTheLoop
number currently of processed method
const Problem * problemPtr
raw pointer to problem
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
MatrixInt facesNodes
nodes on finite element faces
std::string feName
Name of finite element.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
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
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MatrixDouble gaussPts
Matrix of integration points.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
Block data for Nitsche method.
double gamma
Penalty term, see .
string faceElemName
name of element face
double phi
Nitsche method parameter, see .
Range fAces
faces on which constrain is applied
MultiIndexData(int gg, int side, EntityType type)
Common data shared between finite element operators.
std::vector< std::vector< int > > inTetFaceGaussPtsNumber
std::vector< boost::shared_ptr< const NumeredEntFiniteElement > > facesFePtr
std::vector< MatrixDouble > faceNormals
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
std::vector< VectorDouble > rAy
std::vector< MatrixDouble > coordsAtGaussPts
std::map< EntityType, std::vector< std::vector< MatrixDouble > > > dP
derivative of projection matrix in respect DoFs This is EntityType, face, gauss point at face....
std::vector< std::vector< MatrixDouble > > P
projection matrix
std::vector< VectorDouble > cOords
std::vector< MatrixDouble > hoCoordsAtGaussPts
std::vector< MatrixDouble > faceGaussPts
std::vector< EntityHandle > fAces
int getRule(int order)
MyFace(MoFEM::Interface &m_field)
Definition of volume element.
MyVolumeFE(MoFEM::Interface &m_field, BlockData &block_data, CommonData &common_data)
virtual MoFEMErrorCode doAdditionalJobWhenGuassPtsAreCalulated()
MoFEMErrorCode setGaussPts(int order)
Basic operated shared between all Nitsche operators.
virtual MoFEMErrorCode getFaceRadius(int ff)
OpBasicCommon(const std::string field_name, BlockData &nitsche_block_data, CommonData &nitsche_common_data, bool field_disp, const char type)
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)
MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &data, int gg, MatrixDouble &jac)
virtual MoFEMErrorCode calculateP(int gg, int fgg, int ff)
NonlinearElasticElement::BlockData & dAta
MoFEMErrorCode getTractionVariance(int gg, int fgg, int ff, MatrixDouble &jac, MatrixDouble &trac)
Get integration pts data on face.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetFaceData(CommonData &common_data)
Calculate Nitsche method terms on left hand side.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
std::vector< MatrixDouble > kMatrixFace
std::vector< MatrixDouble > kMatrixFace0
OpLhsNormal(const std::string field_name, BlockData &nitsche_block_data, CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, bool field_disp)
std::vector< MatrixDouble > kMatrixFace1
Calculate Nitsche method terms on right hand side.
OpRhsNormal(const std::string field_name, BlockData &nitsche_block_data, CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, bool field_disp)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
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
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< MatrixDouble3by3 > sTress
void tricircumcenter3d_tp(a, b, c, circumcenter, double *xi, double *eta)