OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
prescribedgradientbcweak.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
36 #include "classfactory.h"
37 #include "node.h"
38 #include "masterdof.h"
39 #include "element.h"
40 #include "feinterpol.h"
41 #include "feinterpol2d.h"
42 #include "gausspoint.h"
43 #include "sparsemtrx.h"
46 #include "spatiallocalizer.h"
47 #include "geometry.h"
48 #include "dynamicinputrecord.h"
49 #include "timestep.h"
50 #include "function.h"
51 #include "engngm.h"
52 #include "mathfem.h"
53 #include "sparselinsystemnm.h"
54 #include "unknownnumberingscheme.h"
55 #include "../sm/Materials/structuralmaterial.h"
56 #include "../sm/EngineeringModels/staticstructural.h"
57 
58 #include "timer.h"
59 
60 #include "xfem/XFEMDebugTools.h"
61 
62 #include <cmath>
63 #include <string>
64 #include <sstream>
65 #include <iterator>
66 
67 namespace oofem {
68 
70 {
71  rows.clear();
72 
73  IntArray tracElRows, trac_loc_r;
74 
75  mFirstNode->giveLocationArray({Trac_u, Trac_v}, trac_loc_r, s);
76  tracElRows.followedBy(trac_loc_r);
77  rows = tracElRows;
78 }
79 
81 {
82 
83  Element *dummyEl = NULL;
85 
86  // Material mode should not matter here.
87  MaterialMode matMode = _PlaneStrain;
88 
89  int numPointsPerSeg = 1;
90  mIntRule->SetUpPointsOnLine(numPointsPerSeg, matMode);
91 
92 
93 }
94 
95 
99  mTractionDofIDs( {Trac_u, Trac_v} ),
100  mDispLockDofIDs( {LMP_u, LMP_v} ),
101  mRegularDispDofIDs( {D_u, D_v} ),
105  mMeshIsPeriodic(false),
106  mDuplicateCornerNodes(false),
107  mTangDistPadding(0.0),
108  mTracDofScaling(1.0e0),
109  mpDisplacementLock(NULL),
110  mLockNodeInd(0),
111  mDispLockScaling(1.0),
112  mSpringNodeInd1(-1),
113  mSpringNodeInd2(-1),
114  mSpringPenaltyStiffness(1.0e-3),
115  mPeriodicityNormal({0.0, 1.0}),
116  mDomainSize(0.0),
117  mMirrorFunction(0)
118 {
119  if ( d ) {
120  // Compute bounding box of the domain
122  }
123 
124 }
125 
127 {
128  clear();
129 }
130 
132 {
133 
134  for ( size_t i = 0; i < mpTracElNew.size(); i++ ) {
135  if ( mpTracElNew [ i ] != NULL ) {
136  delete mpTracElNew [ i ];
137  mpTracElNew [ i ] = NULL;
138  }
139  }
140  mpTracElNew.clear();
141 
142  if ( mpDisplacementLock != NULL ) {
143  delete mpDisplacementLock;
144  mpDisplacementLock = NULL;
145  }
146 
147 }
148 
149 //#define DAMAGE_TEST
150 
152 {
153  int numDMan = mpTracElNew.size();
154 
155  if ( mpDisplacementLock != NULL ) {
156  numDMan++;
157  }
158 
159  return numDMan;
160 }
161 
163 {
164  if ( i - 1 < int( mpTracElNew.size() ) ) {
165  return mpTracElNew [ i - 1 ]->mFirstNode.get();
166  } else {
167  OOFEM_ERROR("return mpDisplacementLock")
168  return mpDisplacementLock;
169  }
170 }
171 
173 {
174  IRResultType result;
176  if ( result != IRRT_OK ) {
177  return result;
178  }
180  if ( result != IRRT_OK ) {
181  return result;
182  }
183 
185 // printf("mTractionInterpOrder: %d\n", mTractionInterpOrder);
186 
188 // printf("mNumTractionNodesAtIntersections: %d\n", mNumTractionNodesAtIntersections);
189 
191  OOFEM_ERROR("mNumTractionNodesAtIntersections > 1 is not allowed if mTractionInterpOrder == 0.")
192  }
193 
195 // printf("mTractionNodeSpacing: %d\n", mTractionNodeSpacing);
196 
197  int duplicateCorners = 0;
199 // printf("duplicateCorners: %d\n", duplicateCorners);
200 
201  if ( duplicateCorners == 1 ) {
202  mDuplicateCornerNodes = true;
203  } else {
204  mDuplicateCornerNodes = false;
205  }
206 
207  mTangDistPadding = 0.0;
209 // printf("mTangDistPadding: %e\n", mTangDistPadding);
210 
212 // printf("mTracDofScaling: %e\n", mTracDofScaling );
213 
216 // printf("mPeriodicityNormal: "); mPeriodicityNormal.printYourself();
217 
218 
220 // printf("mMirrorFunction: %d\n", mMirrorFunction );
221 
222  if ( mMirrorFunction == 0 ) {
223  mPeriodicityNormal = {0.0, 1.0};
224  }
225 
226  return IRRT_OK;
227 }
228 
230 {
233 
237 
238  if ( mDuplicateCornerNodes ) {
240  } else {
242  }
243 
246 
247  if ( mMirrorFunction > 0 ) {
250  }
251 }
252 
254 {}
255 
257  CharType type, ValueModeType mode,
258  const UnknownNumberingScheme &s, FloatArray *eNorm)
259 {
261 
262  if ( type == ExternalForcesVector ) {
263  // The external force vector is given by
264  // f_ext = int N^trac H . (x - x_c)
265 
266 
267  for ( TracSegArray* el : mpTracElNew ) {
268 
269  FloatArray contrib;
270  computeExtForceElContrib(contrib, *el, dim, tStep);
271 
272  IntArray rows;
273  el->giveTractionLocationArray(rows, type, s);
274  answer.assemble(contrib, rows);
275 
276  }
277 
278  } else if ( type == InternalForcesVector ) {
279 
280  for ( TracSegArray* el : mpTracElNew ) {
281 
282  for ( GaussPoint *gp: *(el->mIntRule.get()) ) {
283 
284  // Contribution on gamma_plus
285  FloatArray contrib_disp, contrib_trac;
286  IntArray disp_loc_array, trac_loc_array;
287  computeIntForceGPContrib(contrib_disp, disp_loc_array, contrib_trac, trac_loc_array, *el, *gp, dim, tStep, gp->giveGlobalCoordinates(), 1.0, mode, type, s);
288  answer.assemble(contrib_disp, disp_loc_array);
289  answer.assemble(contrib_trac, trac_loc_array);
290 
291 
292  // Contribution on gamma_minus
293  contrib_disp.clear(); contrib_trac.clear();
294  disp_loc_array.clear(); trac_loc_array.clear();
295  FloatArray xMinus;
296  this->giveMirroredPointOnGammaMinus(xMinus, gp->giveGlobalCoordinates());
297  computeIntForceGPContrib(contrib_disp, disp_loc_array, contrib_trac, trac_loc_array, *el, *gp, dim, tStep, xMinus, -1.0, mode, type, s);
298  answer.assemble(contrib_disp, disp_loc_array);
299  answer.assemble(contrib_trac, trac_loc_array);
300 
301  }
302  }
303 
304 
305 
306  if ( mpDisplacementLock != NULL ) {
307  IntArray dispLockRows;
309 
310  FloatArray fe_dispLock;
311 
312  int lockNodePlaceInArray = domain->giveDofManPlaceInArray(mLockNodeInd);
313  FloatArray nodeUnknowns;
314  domain->giveDofManager(lockNodePlaceInArray)->giveUnknownVector(nodeUnknowns,this->giveRegularDispDofIDs(), mode, tStep);
315 
316  for ( int i = 0; i < domain->giveNumberOfSpatialDimensions(); i++ ) {
317  fe_dispLock.push_back(nodeUnknowns [ i ]);
318  }
319 
320  fe_dispLock.times(mDispLockScaling);
321 
322  answer.assemble(fe_dispLock, dispLockRows);
323  }
324 
325  }
326 }
327 
329 {
330 
331  oContrib.clear();
332  FloatArray contrib_gp;
333 
334  for ( GaussPoint *gp: *(iEl.mIntRule.get()) ) {
335 
336  // Fetch global coordinate x
337  const FloatArray &x = gp->giveGlobalCoordinates();
338 
339  // Compute H.[x]
340  FloatArray temp;
341  giveBoundaryCoordVector(temp, x);
342 
343  FloatArray Hx;
344  FloatMatrix grad2D = mGradient;
345  grad2D.resizeWithData(2,2);
346  Hx.beProductOf(grad2D, temp);
347 
348 
349  // For now, assume piecewise constant approx
350  FloatArray Ntrac = FloatArray { 1.0*mTracDofScaling };
351 
352  // N-matrix
353  FloatMatrix Nmat;
354  Nmat.beNMatrixOf(Ntrac, iDim);
355 
356 
357  // Assemble contribution to the global vector directly
358  contrib_gp.beTProductOf(Nmat, Hx);
359  double detJ = 0.5 * iEl.giveLength();
360  double loadLevel = this->giveTimeFunction()->evaluateAtTime(tStep->giveTargetTime());
361  contrib_gp.times(-detJ * gp->giveWeight()*loadLevel);
362 
363  oContrib.add(contrib_gp);
364  }
365 
366 }
367 
368 void PrescribedGradientBCWeak :: computeIntForceGPContrib(FloatArray &oContrib_disp, IntArray &oDisp_loc_array, FloatArray &oContrib_trac, IntArray &oTrac_loc_array,TracSegArray &iEl, GaussPoint &iGP, int iDim, TimeStep *tStep, const FloatArray &iBndCoord, const double &iScaleFac, ValueModeType mode, CharType type, const UnknownNumberingScheme &s)
369 {
370 
372 
373  FloatMatrix contrib;
374  assembleTangentGPContributionNew(contrib, iEl, iGP, iScaleFac, iBndCoord);
375 
376  // Compute vector of traction unknowns
377  FloatArray tracUnknowns;
378  iEl.mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), mode, tStep);
379 
380  iEl.giveTractionLocationArray(oTrac_loc_array, type, s);
381 
382  FloatArray dispElLocCoord, closestPoint;
383  Element *dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, iBndCoord );
384 
385  // Compute vector of displacement unknowns
386  FloatArray dispUnknowns;
387  int numDMan = dispEl->giveNumberOfDofManagers();
388  for ( int i = 1; i <= numDMan; i++ ) {
389  FloatArray nodeUnknowns;
390  DofManager *dMan = dispEl->giveDofManager(i);
391 
392  IntArray dispIDs = giveRegularDispDofIDs();
393  if ( domain->hasXfemManager() ) {
395  dispIDs.followedBy(xMan->giveEnrichedDofIDs(*dMan));
396  }
397 
398  dMan->giveUnknownVector(nodeUnknowns, dispIDs,mode, tStep);
399  dispUnknowns.append(nodeUnknowns);
400 
401  }
402 
403  dispEl->giveLocationArray(oDisp_loc_array, s);
404 
405 
406  oContrib_disp.beTProductOf(contrib, tracUnknowns);
407  oContrib_disp.negated();
408 
409  oContrib_trac.beProductOf(contrib, dispUnknowns);
410  oContrib_trac.negated();
411 }
412 
414  CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s, double scale)
415 {
416  std::vector<FloatArray> gpCoordArray;
417 
418  if ( type == TangentStiffnessMatrix || type == SecantStiffnessMatrix || type == ElasticStiffnessMatrix ) {
419 
420  for ( TracSegArray* el : mpTracElNew ) {
421 
422  for ( GaussPoint *gp: *(el->mIntRule.get()) ) {
423 
424  gpCoordArray.push_back( gp->giveGlobalCoordinates() );
425  assembleGPContrib(answer, tStep, type, r_s, c_s, *el, *gp, scale);
426  }
427  }
428 
429  if ( mpDisplacementLock != NULL ) {
431  FloatMatrix KeDispLock(nsd, nsd);
432  KeDispLock.beUnitMatrix();
433  KeDispLock.times(mDispLockScaling);
434 
435  int placeInArray = domain->giveDofManPlaceInArray(mLockNodeInd);
436  DofManager *node = domain->giveDofManager(placeInArray);
437 
438  IntArray lockRows, lockCols, nodeRows, nodeCols;
440  node->giveLocationArray(giveRegularDispDofIDs(), nodeRows, r_s);
441 
443  node->giveLocationArray(giveRegularDispDofIDs(), nodeCols, c_s);
444 
445  answer.assemble(lockRows, nodeCols, KeDispLock);
446  answer.assemble(nodeRows, lockCols, KeDispLock);
447 
448  FloatMatrix KZero( lockRows.giveSize(), lockCols.giveSize() );
449  KZero.zero();
450  KZero.times(scale);
451  answer.assemble(lockRows, lockCols, KZero);
452  }
453 
454 
456  FloatMatrix KeDispLock(nsd, nsd);
457  KeDispLock.beUnitMatrix();
458  KeDispLock.times(mSpringPenaltyStiffness);
459 
460 // printf("mSpringNodeInd1: %d\n", mSpringNodeInd1);
461  int placeInArray = domain->giveDofManPlaceInArray(mSpringNodeInd1);
462  DofManager *node1 = domain->giveDofManager(placeInArray);
463 
464  IntArray nodeRows, nodeCols;
465  node1->giveLocationArray(giveRegularDispDofIDs(), nodeRows, r_s);
466  node1->giveLocationArray(giveRegularDispDofIDs(), nodeCols, c_s);
467  KeDispLock.times(scale);
468  answer.assemble(nodeRows, nodeCols, KeDispLock);
469 
470 
471  } else {
472  printf("Skipping assembly in PrescribedGradientBCWeak::assemble().\n");
473  }
474 
475 // std :: string fileName("TracGpCoord.vtk");
476 // XFEMDebugTools :: WritePointsToVTK(fileName, gpCoordArray);
477 }
478 
480  CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)
481 
482 {
483 // printf("Entering PrescribedGradientBCWeak :: assembleExtraDisplock.\n");
484 
485 #if 1
486 
488  FloatMatrix KeDispLock(nsd, nsd);
489  KeDispLock.zero();
490 
491  // Lock in y-direction to get rid of rigid body rotation.
492  double scaling = 1.0e0;
493  KeDispLock.at(2,2) = mSpringPenaltyStiffness*scaling;
494 
495  IntArray nodeRows, nodeCols;
496 
497 // printf("mSpringNodeInd2: %d\n", mSpringNodeInd2);
498  int placeInArray = domain->giveDofManPlaceInArray(mSpringNodeInd2);
499  DofManager *node1 = domain->giveDofManager(placeInArray);
500 
501  node1->giveLocationArray(giveRegularDispDofIDs(), nodeRows, r_s);
502  node1->giveLocationArray(giveRegularDispDofIDs(), nodeCols, c_s);
503 
504  answer.assemble(nodeRows, nodeCols, KeDispLock);
505 
506 
507 #else
509  FloatMatrix KeDispLock(2*nsd, 2*nsd);
510  KeDispLock.zero();
511 
512  // Lock in y-direction to get rid of rigid body rotation.
513  double scaling = 1.0e0;
514  KeDispLock.at(2,2) = mSpringPenaltyStiffness*scaling;
515  KeDispLock.at(2,nsd+1) = -mSpringPenaltyStiffness*scaling;
516  KeDispLock.at(nsd+1,2) = -mSpringPenaltyStiffness*scaling;
517  KeDispLock.at(nsd+1,nsd+1) = mSpringPenaltyStiffness*scaling;
518 // KeDispLock.times(mSpringPenaltyStiffness);
519 
520  IntArray nodeRows, nodeCols;
521 
522 
523  int placeInArray = domain->giveDofManPlaceInArray(mSpringNodeInd2);
524  DofManager *node1 = domain->giveDofManager(placeInArray);
525 
526  IntArray nodeRows1, nodeCols1;
527  node1->giveLocationArray(giveRegularDispDofIDs(), nodeRows1, r_s);
528  node1->giveLocationArray(giveRegularDispDofIDs(), nodeCols1, c_s);
529 
530  nodeRows.followedBy(nodeRows1);
531  nodeCols.followedBy(nodeCols1);
532 
533 
535  DofManager *node2 = domain->giveDofManager(placeInArray);
536 
537  IntArray nodeRows2, nodeCols2;
538  node2->giveLocationArray(giveRegularDispDofIDs(), nodeRows2, r_s);
539  node2->giveLocationArray(giveRegularDispDofIDs(), nodeCols2, c_s);
540 
541  nodeRows.followedBy(nodeRows2);
542  nodeCols.followedBy(nodeCols2);
543 
544  answer.assemble(nodeRows, nodeCols, KeDispLock);
545 #endif
546 
547 }
548 
550  CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s, TracSegArray &iEl, GaussPoint &iGP, double k)
551 {
552 
554 
556  // Gamma_plus
557  FloatMatrix contrib;
558  assembleTangentGPContributionNew(contrib, iEl, iGP, -1.0, iGP.giveGlobalCoordinates());
559 
560  // Compute vector of traction unknowns
561  FloatArray tracUnknowns;
562  iEl.mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);
563 
564  IntArray trac_rows;
565  iEl.giveTractionLocationArray(trac_rows, type, r_s);
566 
567 
568  FloatArray dispElLocCoord, closestPoint;
569  Element *dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, iGP.giveGlobalCoordinates() );
570 
571  IntArray disp_cols;
572  dispEl->giveLocationArray(disp_cols, c_s);
573 
574  contrib.times(k);
575  answer.assemble(trac_rows, disp_cols, contrib);
576 
577  FloatMatrix contribT;
578  contribT.beTranspositionOf(contrib);
579  answer.assemble(disp_cols, trac_rows, contribT);
580 
581 
582 
584  // Gamma_minus
585  contrib.clear();
586  FloatArray xMinus;
588  assembleTangentGPContributionNew(contrib, iEl, iGP, 1.0, xMinus);
589 
590  // Compute vector of traction unknowns
591  tracUnknowns.clear();
592  iEl.mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);
593 
594  trac_rows.clear();
595  iEl.giveTractionLocationArray(trac_rows, type, r_s);
596 
597 
598  dispElLocCoord.clear(); closestPoint.clear();
599  dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, xMinus );
600 
601  disp_cols.clear();
602  dispEl->giveLocationArray(disp_cols, c_s);
603  contrib.times(k);
604  answer.assemble(trac_rows, disp_cols, contrib);
605 
606  contribT.clear();
607  contribT.beTranspositionOf(contrib);
608  answer.assemble(disp_cols, trac_rows, contribT);
609 
610 
611  // Assemble zeros on diagonal (required by PETSc solver)
612  FloatMatrix KZero(1,1);
613  KZero.zero();
614  for ( int i : trac_rows) {
615  answer.assemble(IntArray({i}), IntArray({i}), KZero);
616  }
617 }
618 
619 void PrescribedGradientBCWeak :: giveLocationArrays(std :: vector< IntArray > &rows, std :: vector< IntArray > &cols, CharType type,
620  const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)
621 {
622 
623 }
624 
626  const UnknownNumberingScheme &s)
627 {
628  // Used for the condensation when computing the macroscopic tangent.
629 
630 
631  rows.clear();
632 
633  // Loop over traction elements
634  for ( TracSegArray* el : mpTracElNew ) {
635 
636  IntArray trac_loc_r;
637  el->mFirstNode->giveLocationArray(giveTracDofIDs(), trac_loc_r, s);
638 
639  rows.followedBy(trac_loc_r);
640  }
641 
642 #if 0
643  rows.clear();
644 
645  // Loop over traction elements
646  for ( size_t tracElInd = 0; tracElInd < mpTractionElements.size(); tracElInd++ ) {
647  IntArray tracElRows, trac_loc_r;
648 
649  const TractionElement &tEl = * ( mpTractionElements [ tracElInd ] );
650  for ( int tracNodeInd : tEl.mTractionNodeInd ) {
651  Node *tNode = mpTractionNodes [ tracNodeInd ];
652  tNode->giveLocationArray(giveTracDofIDs(), trac_loc_r, s);
653  tracElRows.followedBy(trac_loc_r);
654  }
655 
656  rows.followedBy(tracElRows);
657  }
658 
659 
660  if ( mpDisplacementLock != NULL ) {
661  IntArray dispLock_r;
663 
664  rows.followedBy(dispLock_r);
665  }
666 #endif
667 }
668 
670  const UnknownNumberingScheme &s)
671 {
672  // Used for the condensation when computing the macroscopic tangent.
673 
674 }
675 
677 {
678  //double rve_size = this->domainSize();
679  const int dim = domain->giveNumberOfSpatialDimensions();
680 
681  IntArray loc_t;
683  giveTractionLocationArray(loc_t, fnum);
684 
685  int num_t_eq = loc_t.giveSize();
686 
687 // o_x_times_N.resize(3, num_t_eq);
688  o_x_times_N.resize(num_t_eq,3);
689 
690  IntArray cols = {1,2,3};
691 
692  int trac_el_ind = 1;
693 
694  for ( auto *el: mpTracElNew ) {
695 
696  IntArray rows = {trac_el_ind*2-1,trac_el_ind*2};
697 
698  for ( GaussPoint *gp: *(el->mIntRule.get()) ) {
699 
700  FloatMatrix contrib(2,3);
701 
702  // For now, assume piecewise constant approx
703  FloatArray Ntrac = FloatArray { 1.0*mTracDofScaling };
704 
705  // N-matrix
706  FloatMatrix Nmat;
707  Nmat.beNMatrixOf(Ntrac, dim);
708 // printf("Nmat: "); Nmat.printYourself();
709 
710  // Fetch global coordinate x
711  const FloatArray &x = gp->giveGlobalCoordinates();
712 
713 // // Compute vector of traction unknowns
714 // FloatArray tracUnknowns;
715 // el->mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);
716 
717 
718 // FloatArray traction;
719 // traction.beProductOf(Nmat, tracUnknowns);
720 
721  FloatArray tmp;
722  giveBoundaryCoordVector(tmp, x);
723 
724  FloatMatrix coord_mat(2,3);
725  coord_mat.at(1,1) = tmp.at(1);
726  coord_mat.at(2,2) = tmp.at(2);
727 
728 // coord_mat.at(3,2) = tmp.at(1);
729 
730 // coord_mat.at(3,1) = tmp.at(2);
731 // coord_mat.at(4,2) = tmp.at(1);
732 
733 // coord_mat.at(4,1) = tmp.at(2);
734 // coord_mat.at(3,2) = tmp.at(1);
735 
736  coord_mat.at(1,3) = 0.5*tmp.at(2);
737  coord_mat.at(2,3) = 0.5*tmp.at(1);
738 
739 // coord_mat.at(4,2) = 0.5*tmp.at(2);
740 // coord_mat.at(4,1) = 0.5*tmp.at(1);
741 
742 
743 // FloatMatrix contrib;
744 // contrib.beDyadicProductOf(traction, tmp);
745 
746  contrib.beProductOf(Nmat, coord_mat);
747 
748 // contrib = coord_mat;
749 
750  double detJ = 0.5 * el->giveLength();
751  contrib.times( detJ * gp->giveWeight() );
752 
753 // printf("\n\ncontrib: "); contrib.printYourself();
754 // printf("rows: "); rows.printYourself();
755 // printf("cols: "); cols.printYourself();
756 
757  o_x_times_N.assemble(contrib, rows, cols);
758  }
759 
760  trac_el_ind++;
761  }
762 }
763 
765 {
766  //double rve_size = this->domainSize();
767  const int dim = domain->giveNumberOfSpatialDimensions();
768 
769  IntArray loc_t;
771  giveTractionLocationArray(loc_t, fnum);
772 
773  int num_t_eq = loc_t.giveSize();
774 
775 // o_x_times_N.resize(4, num_t_eq);
776  o_x_times_N.resize(3, num_t_eq);
777 
778  IntArray rows = {1,2,3};
779 
780  int trac_el_ind = 1;
781 
782  for ( auto *el: mpTracElNew ) {
783 
784  IntArray cols = {trac_el_ind*2-1,trac_el_ind*2};
785 
786  for ( GaussPoint *gp: *(el->mIntRule.get()) ) {
787 
788  FloatMatrix contrib(4,2);
789 
790  // For now, assume piecewise constant approx
791  FloatArray Ntrac = FloatArray { 1.0*mTracDofScaling };
792 
793  // N-matrix
794  FloatMatrix Nmat;
795  Nmat.beNMatrixOf(Ntrac, dim);
796 // printf("Nmat: "); Nmat.printYourself();
797 
798  // Fetch global coordinate x
799  const FloatArray &x = gp->giveGlobalCoordinates();
800 
801 // // Compute vector of traction unknowns
802 // FloatArray tracUnknowns;
803 // el->mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);
804 
805 
806 // FloatArray traction;
807 // traction.beProductOf(Nmat, tracUnknowns);
808 
809  FloatArray tmp;
810  giveBoundaryCoordVector(tmp, x);
811 
812  FloatMatrix coord_mat(4,2);
813  coord_mat.at(1,1) = tmp.at(1);
814  coord_mat.at(2,2) = tmp.at(2);
815 
816 // coord_mat.at(3,2) = tmp.at(1);
817 // coord_mat.at(3,1) = tmp.at(2);
818 
819 // coord_mat.at(3,1) = tmp.at(2);
820 // coord_mat.at(4,2) = tmp.at(1);
821 
822  coord_mat.at(3,1) = 0.5*tmp.at(2);
823  coord_mat.at(3,2) = 0.5*tmp.at(1);
824 //
825 // coord_mat.at(4,2) = 0.5*tmp.at(2);
826 // coord_mat.at(4,1) = 0.5*tmp.at(1);
827 
828 
829 // FloatMatrix contrib;
830 // contrib.beDyadicProductOf(traction, tmp);
831 
832  contrib.beProductOf(coord_mat, Nmat);
833 
834 // contrib = coord_mat;
835 
836  double detJ = 0.5 * el->giveLength();
837  contrib.times( detJ * gp->giveWeight() );
838 
839 // printf("\n\ncontrib: "); contrib.printYourself();
840 // printf("rows: "); rows.printYourself();
841 // printf("cols: "); cols.printYourself();
842 
843  o_x_times_N.assemble(contrib, rows, cols);
844  }
845 
846  trac_el_ind++;
847  }
848 }
849 
850 
852 {
853  double Lx = mUC[0] - mLC[0];
854  double Ly = mUC[1] - mLC[1];
855  double dSize = Lx*Ly;
856 // printf("dSize: %e\n", dSize);
857 
858  const int dim = domain->giveNumberOfSpatialDimensions();
859  FloatMatrix stressMatrix(dim, dim);
860 
861  for ( auto *el: mpTracElNew ) {
862 
863  for ( GaussPoint *gp: *(el->mIntRule.get()) ) {
864 
865  // For now, assume piecewise constant approx
866  FloatArray Ntrac = FloatArray { 1.0*mTracDofScaling };
867 
868  // N-matrix
869  FloatMatrix Nmat;
870  Nmat.beNMatrixOf(Ntrac, dim);
871 
872  // Fetch global coordinate x
873  const FloatArray &x = gp->giveGlobalCoordinates();
874 
875  // Compute vector of traction unknowns
876  FloatArray tracUnknowns;
877  el->mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);
878 
879 
880  FloatArray traction;
881  traction.beProductOf(Nmat, tracUnknowns);
882 
883  FloatArray tmp;
884  giveBoundaryCoordVector(tmp, x);
885 
886  FloatMatrix contrib;
887  contrib.beDyadicProductOf(traction, tmp);
888 
889  double detJ = 0.5 * el->giveLength();
890  contrib.times( detJ * gp->giveWeight() );
891 
892  for ( int m = 0; m < dim; m++ ) {
893  for ( int n = 0; n < dim; n++ ) {
894  stressMatrix(m, n) += contrib(m, n);
895  }
896  }
897 
898  }
899 
900  }
901 
902  if ( dim == 2 ) {
903  sigma = {
904  stressMatrix(0, 0), stressMatrix(1, 1), 0.0, 0.0, 0.0, 0.5*(stressMatrix(0, 1) + stressMatrix(1, 0))
905  };
906  } else {
907  sigma.beVectorForm(stressMatrix);
908  }
909 
910  sigma.times(1.0 / dSize);
911 
912 }
913 
914 //#define TIME_INFO
915 
917 {
918 #ifdef TIME_INFO
919  static double tot_time = 0.0;
920  Timer timer;
921  timer.startTimer();
922 
923  static double assemble_time = 0.0;
924  Timer assemble_timer;
925 
926 #endif
927 
928  E.resize(9,9);
929 
930 
931  // Extract the relevant submatrices from the RVE problem.
932  // At equilibrium, the RVE problem has the structure
933  //
934  // [ S C; C^T 0 ]*[a_u a_t] = [0; f]
935  //
936  // where a_u and a_t denote displacement and traction dofs, respectively.
937  // We need to extract S and C.
938 
939  EngngModel *rve = this->giveDomain()->giveEngngModel();
941  std :: unique_ptr< SparseLinearSystemNM > solver(
942  classFactory.createSparseLinSolver( ST_Petsc, this->domain, this->domain->giveEngngModel() ) ); // = rve->giveLinearSolver();
943  bool symmetric_matrix = false;
944  SparseMtrxType stype = solver->giveRecommendedMatrix(symmetric_matrix);
945 // double rve_size = this->domainSize();
946  double Lx = mUC[0] - mLC[0];
947  double Ly = mUC[1] - mLC[1];
948  double rve_size = Lx*Ly;
949 
950 
952  std :: unique_ptr< SparseMtrx > Kmicro( classFactory.createSparseMtrx(stype) );
953  if ( !Kmicro ) {
954  OOFEM_ERROR("Couldn't create sparse matrix of type %d\n", stype);
955  }
956 
957  StaticStructural *rveStatStruct = dynamic_cast<StaticStructural*>(rve);
958  if ( rveStatStruct ) {
959 // printf("Successfully casted rve to StaticStructural.\n");
960 
961  if ( rveStatStruct->stiffnessMatrix ) {
962  Kmicro.reset( rveStatStruct->stiffnessMatrix->GiveCopy() );
963  }
964  }
965 
966 
967 #ifdef TIME_INFO
968  assemble_timer.startTimer();
969 #endif
970 
971  if ( Kmicro->giveNumberOfColumns() == 0 ) {
972 // printf("Rebuilding stiffness matrix.\n");
973  Kmicro->buildInternalStructure(rve, this->domain->giveNumber(), fnum);
974  rve->assemble(*Kmicro, tStep, TangentAssembler(TangentStiffness), fnum, this->domain);
975  }
976 // else {
977 // printf("Using existing stiffness matrix.\n");
978 // }
979 
980  assembleExtraDisplock(*Kmicro, tStep, TangentStiffnessMatrix, fnum, fnum);
981 #ifdef TIME_INFO
982  assemble_timer.stopTimer();
983  assemble_time += assemble_timer.getUtime();
984  printf("Assembly time for RVE tangent: %e\n", assemble_time);
985 #endif
986 
987 
988  // Fetch displacement and traction location arrays
989  IntArray loc_u, loc_t;
990 
991  giveTractionLocationArray(loc_t, fnum);
992 
993  int neq = Kmicro->giveNumberOfRows();
994  loc_u.resize(neq - loc_t.giveSize());
995  int k = 1;
996  for ( int i = 1; i <= neq; i++ ) {
997  if ( !loc_t.contains(i) ) {
998  loc_u.at(k) = i;
999  k++;
1000  }
1001  }
1002 
1003 
1004  // Fetch the submatrices
1005  std :: unique_ptr< SparseMtrx > S(Kmicro->giveSubMatrix(loc_u, loc_u));
1006  // NOTE: Kus is actually a dense matrix, but we have to make it a dense matrix first
1007  std :: unique_ptr< SparseMtrx > C(Kmicro->giveSubMatrix(loc_u, loc_t));
1008  FloatMatrix Cd;
1009  C->toFloatMatrix(Cd);
1010 
1011 
1012  // Let Sm = C. Solve for m.
1013  FloatMatrix m;
1014  solver->solve(*S, Cd, m);
1015 
1016 
1017  // Compute G := C^T m. (Which could, formally, be written as G = C^T S^-1 C.)
1018  FloatMatrix G;
1019  G.beTProductOf(Cd, m);
1020 
1021 
1022  // Compute D := \int x \otimes N
1023  FloatMatrix D;
1025 // FloatMatrix D;
1026 // compute_x_times_N_1(DT2);
1027 // D.beTranspositionOf(DT);
1028 
1029 
1030  // Let Gp = D. Solve for p. (Which could, formally, be written as p = G^-1 D = (C^T S^-1 C)^-1 D.)
1031  // Need to make G sparse;
1032  std :: unique_ptr< SparseMtrx > Gs( classFactory.createSparseMtrx(stype) );
1033  if ( !Gs ) {
1034  OOFEM_ERROR("Couldn't create sparse matrix of type %d\n", stype);
1035  }
1036 
1037  int num_eq_G = G.giveNumberOfRows();
1038  IntArray loc_G(num_eq_G);
1039 
1040  for ( int i = 1; i <= num_eq_G; i++) {
1041  loc_G.at(i) = i;
1042  }
1043 
1044  Gs->buildInternalStructure(rve, num_eq_G, num_eq_G, loc_G, loc_G);
1045  Gs->assemble(loc_G, loc_G, G);
1046 
1047  Gs->assembleBegin();
1048  Gs->assembleEnd();
1049 
1050 // Gs->writeToFile("Gs.txt");
1051 
1052 
1053  FloatMatrix p;
1054  solver->solve(*Gs, D, p);
1055 
1056 
1057  // Compute d_sigma_depsilon = D^T p.
1058  FloatMatrix Ered;
1059  Ered.beTProductOf(D,p);
1060 // rve_size = 25.0;
1061  Ered.times(1.0/rve_size);
1062 // printf("rve_size: %e\n", rve_size);
1063 // Ered.printYourself();
1064 
1065  IntArray indx;
1066  StructuralMaterial :: giveVoigtVectorMask(indx, _PlaneStress);
1067 
1068 // FloatMatrix EredT;
1069 // EredT.beTranspositionOf(Ered);
1070 // Ered.add(EredT);
1071 // Ered.times(0.5);
1072 
1073 // Ered.at(1,3) = 0.0;
1074 // Ered.at(2,3) = 0.0;
1075 // Ered.at(3,1) = 0.0;
1076 // Ered.at(3,2) = 0.0;
1077 
1078  E.assemble(Ered, indx, indx);
1079 // E.printYourself();
1080 
1081 #ifdef TIME_INFO
1082  timer.stopTimer();
1083  tot_time += timer.getUtime();
1084 // printf("Total time for RVE tangent: %e\n", tot_time);
1085 #endif
1086 
1087 }
1088 
1089 void PrescribedGradientBCWeak :: giveTractionElNormal(size_t iElInd, FloatArray &oNormal, FloatArray &oTangent) const
1090 {
1091  FloatArray xS, xE;
1092  giveTractionElCoord(iElInd, xS, xE);
1093 
1094  oTangent.beDifferenceOf(xE, xS);
1095  oTangent.normalize();
1096 
1097  oNormal = {
1098  oTangent [ 1 ], -oTangent [ 0 ]
1099  };
1100 }
1101 
1102 void PrescribedGradientBCWeak :: giveTractionElArcPos(size_t iElInd, double &oXiStart, double &oXiEnd) const
1103 {
1104  FloatArray xS, xE;
1105  giveTractionElCoord(iElInd, xS, xE);
1106 
1107  FloatArray xC;
1108  xC.beScaled(0.5, xS);
1109  xC.add(0.5, xE);
1110  int sideIndex = giveSideIndex(xC);
1111 
1112  const double nodeDistTol = 1.0e-15;
1113  ArcPosSortFunction3< bool >sortFunc(mLC, mUC, nodeDistTol, sideIndex);
1114 
1115  oXiStart = sortFunc.calcArcPos(xS);
1116  oXiEnd = sortFunc.calcArcPos(xE);
1117 }
1118 
1120 {
1121  Set *setPointer = this->giveDomain()->giveSet(this->set);
1122  oBoundaries = setPointer->giveBoundaryList();
1123 }
1124 
1125 void PrescribedGradientBCWeak :: giveTraction(size_t iElInd, FloatArray &oStartTraction, FloatArray &oEndTraction, ValueModeType mode, TimeStep *tStep)
1126 {
1127  // For now, assuming piecewise constant traction
1128  mpTracElNew[iElInd]->mFirstNode->giveUnknownVector(oStartTraction, giveTracDofIDs(), mode, tStep);
1129  oStartTraction.times(mTracDofScaling);
1130  mpTracElNew[iElInd]->mFirstNode->giveUnknownVector(oEndTraction, giveTracDofIDs(), mode, tStep);
1131  oEndTraction.times(mTracDofScaling);
1132 }
1133 
1135 {
1136 // printf("Recomputing traction mesh.\n");
1137  clear();
1138  postInitialize();
1139 }
1140 
1141 void PrescribedGradientBCWeak :: createTractionMesh(bool iEnforceCornerPeriodicity, int iNumSides)
1142 {
1143  bool split_at_holes = true;
1144 
1145  const double l_s = mUC[0] - mLC[0];
1146  const double minPointDist = 1.0e-4*l_s;
1147 
1148  // Find holes intersecting the RVE boundary so that these can be excluded
1149  std::vector<FloatArray> holeCoordUnsorted, allCoordUnsorted;
1150  findHoleCoord(holeCoordUnsorted, allCoordUnsorted);
1151 
1152 
1153  // Add corner points
1154  holeCoordUnsorted.push_back( {mUC[0], mLC[1]} );
1155  allCoordUnsorted.push_back( {mUC[0], mLC[1]} );
1156 
1157  holeCoordUnsorted.push_back( {mUC[0], mUC[1]} );
1158  allCoordUnsorted.push_back( {mUC[0], mUC[1]} );
1159 
1160  holeCoordUnsorted.push_back( {mLC[0], mUC[1]} );
1161  allCoordUnsorted.push_back( {mLC[0], mUC[1]} );
1162 
1163 
1164  // Add crack-boundary intersections
1165  findCrackBndIntersecCoord(holeCoordUnsorted);
1166 
1167 
1168 
1169  // Add periodicity points
1170  findPeriodicityCoord(holeCoordUnsorted);
1171 
1172 
1173  // Sort arrays in terms of arc length along the RVE boundary
1174  std :: sort( holeCoordUnsorted.begin(), holeCoordUnsorted.end(), ArcPosSortFunction4( mLC, mUC, 1.0e-4 ) );
1175  std :: sort( allCoordUnsorted.begin(), allCoordUnsorted.end(), ArcPosSortFunction4( mLC, mUC, 1.0e-4 ) );
1176 
1177 
1178  // Add refinement points
1179  int pointsPassed = 0;
1180  for ( const auto &x : allCoordUnsorted ) {
1181 
1182  if ( pointsPassed >= mTractionNodeSpacing ) {
1183  holeCoordUnsorted.push_back(x);
1184  pointsPassed = 0;
1185  }
1186 
1187 
1188  pointsPassed++;
1189  }
1190 
1191 
1192  // Sort again
1193  std :: sort( holeCoordUnsorted.begin(), holeCoordUnsorted.end(), ArcPosSortFunction4( mLC, mUC, 1.0e-4 ) );
1194  std :: sort( allCoordUnsorted.begin(), allCoordUnsorted.end(), ArcPosSortFunction4( mLC, mUC, 1.0e-4 ) );
1195 
1196 
1197 
1198  // Remove points that are too close to each other
1199  removeClosePoints(holeCoordUnsorted, minPointDist);
1200 
1201  removeClosePoints(allCoordUnsorted, minPointDist);
1202 
1203 
1204  // Create two arrays of segments, where each array represents the coarsest possible traction
1205  // mesh on one side of the RVE
1206  ArcPosSortFunction4 arcPosFunc( mLC, mUC, 1.0e-4 );
1207 
1208  std :: vector< TracSegArray * > tracElNew0, tracElNew1;
1209  tracElNew0.push_back(new TracSegArray());
1210  tracElNew1.push_back(new TracSegArray());
1211 
1212  for (size_t i = 1; i < holeCoordUnsorted.size(); i++) {
1213 
1214  FloatArray xS = holeCoordUnsorted[i-1];
1215  xS.resizeWithValues(2);
1216  FloatArray xE = holeCoordUnsorted[i];
1217  xE.resizeWithValues(2);
1218  const FloatArray xC = {0.5*(xS[0]+xE[0]), 0.5*(xS[1]+xE[1])};
1219 
1220  if ( arcPosFunc.calcArcPos(xC) < 2.*l_s ) {
1221  tracElNew0[0]->mInteriorSegments.push_back( Line(xS, xE) );
1222  } else {
1223  tracElNew1[0]->mInteriorSegments.push_back( Line(xS, xE) );
1224  }
1225  }
1226 
1227  // Remove segments located in holes
1228  removeSegOverHoles(*(tracElNew0[0]), 1.0e-4);
1229  removeSegOverHoles(*(tracElNew1[0]), 1.0e-4);
1230 
1231  if ( split_at_holes ) {
1232  splitSegments(tracElNew0);
1233  splitSegments(tracElNew1);
1234  }
1235 
1236  // Identify additional points that can be used to refine the traction mesh
1237 
1238 
1239 
1240 
1242  // Create traction dofs
1243  int numNodes = domain->giveNumberOfDofManagers();
1244  int totNodesCreated = 0;
1245 
1246 
1247  // For each side (0 and 1), loop over over elements
1248 
1249  // We may always create the first node on the element.
1250  // For the linear approximation, it may need to be a slave node,
1251  // depending on which element it is.
1252 
1253  // For now, consider only piecewise constant approximations. Then,
1254  // we can always create on node on each element.
1255 
1256  // RVE side at x=L
1257  for ( TracSegArray* el : tracElNew0 ) {
1258 
1260  // Create first node
1261  totNodesCreated++;
1262 
1263  el->mFirstNode.reset( new Node(numNodes + 1, domain) );
1264  el->mFirstNode->setGlobalNumber(numNodes + 1);
1265  for ( auto &dofId: giveTracDofIDs() ) {
1266  el->mFirstNode->appendDof( new MasterDof(el->mFirstNode.get(), ( DofIDItem ) dofId) );
1267  }
1268 
1269  el->mFirstNode->setCoordinates( el->mInteriorSegments[0].giveVertex(1) );
1270  numNodes++;
1271  }
1272 
1273 
1274  // RVE side at y=L
1275  for ( TracSegArray* el : tracElNew1 ) {
1276 
1278  // Create first node
1279  totNodesCreated++;
1280 
1281  el->mFirstNode.reset( new Node(numNodes + 1, domain) );
1282  el->mFirstNode->setGlobalNumber(numNodes + 1);
1283  for ( auto &dofId: giveTracDofIDs() ) {
1284  el->mFirstNode->appendDof( new MasterDof(el->mFirstNode.get(), ( DofIDItem ) dofId) );
1285  }
1286 
1287  el->mFirstNode->setCoordinates( el->mInteriorSegments[0].giveVertex(1) );
1288  numNodes++;
1289  }
1290 
1291  if ( mMeshIsPeriodic && false ) {
1292  // Lock displacement in one node if we use periodic BCs
1293 
1294  int numNodes = domain->giveNumberOfDofManagers();
1295  mpDisplacementLock = new Node(numNodes + 1, domain);
1297 
1298 
1299  for ( auto &dofid: giveDispLockDofIDs() ) {
1301  }
1302  }
1303 
1304 
1305  // Nodes to lock in order to prevent rigid body motion
1307  FloatArray x1 = mLC;
1308 // printf("x1: "); x1.printYourself();
1309  double maxDist = 1.0e10;
1310  Node *node1 = localizer->giveNodeClosestToPoint(x1, maxDist);
1311  mSpringNodeInd1 = node1->giveGlobalNumber();
1312 // printf("mSpringNodeInd1: %d\n", mSpringNodeInd1 );
1313 
1314  FloatArray x2 = {mUC.at(1), mLC.at(2)};
1315 // FloatArray x2 = {mUC.at(1), mUC.at(2)};
1316 // printf("x2: "); x2.printYourself();
1317  Node *node2 = localizer->giveNodeClosestToPoint(x2, maxDist);
1318  mSpringNodeInd2 = node2->giveGlobalNumber();
1319 // printf("mSpringNodeInd2: %d\n", mSpringNodeInd2 );
1320 
1321 
1322  FloatArray x3 = {mLC.at(1), mUC.at(2)};
1323 // FloatArray x2 = {mUC.at(1), mUC.at(2)};
1324 // printf("x3: "); x3.printYourself();
1325  Node *node3 = localizer->giveNodeClosestToPoint(x3, maxDist);
1326  mSpringNodeInd3 = node3->giveGlobalNumber();
1327 // printf("mSpringNodeInd3: %d\n", mSpringNodeInd3 );
1328 
1329 
1330  mpTracElNew.reserve(tracElNew0.size() + tracElNew1.size());
1331  std::move(tracElNew0.begin(), tracElNew0.end(), std::inserter(mpTracElNew, mpTracElNew.end()));
1332  std::move(tracElNew1.begin(), tracElNew1.end(), std::inserter(mpTracElNew, mpTracElNew.end()));
1333  tracElNew0.clear();
1334  tracElNew1.clear();
1335 
1336 
1338  // Segment arrays for Gauss quadrature
1339  size_t i = 0;
1340 
1341  for ( TracSegArray* el : mpTracElNew ) {
1342 
1343  const FloatArray &xS = el->mInteriorSegments[0].giveVertex(1);
1344  const double arcPosXS = arcPosFunc.calcArcPos(xS);
1345 
1346  const FloatArray &xE = el->mInteriorSegments.back().giveVertex(2);
1347  const double arcPosXE = arcPosFunc.calcArcPos(xE);
1348 
1349  while (i < allCoordUnsorted.size()) {
1350 
1351  FloatArray x = allCoordUnsorted[i];
1352  x.resizeWithValues(2);
1353  const double arcPosX = arcPosFunc.calcArcPos(x);
1354 
1355  if ( arcPosX > (arcPosXS+minPointDist) && arcPosX < (arcPosXE-minPointDist) ) {
1356  el->mInteriorSegmentsPointsFine.push_back(x);
1357  }
1358 
1359  if ( arcPosX > arcPosXE ) {
1360  break;
1361  }
1362 
1363  i++;
1364  }
1365  }
1366 
1367  // Now we have the necessary points on each traction element.
1368  // The next step is to create splitted segments.
1369  for ( TracSegArray* el : mpTracElNew ) {
1370 
1371  i = 0;
1372 
1373  for ( Line &line : el->mInteriorSegments ) {
1374  FloatArray xS = line.giveVertex(1);
1375  xS.resizeWithValues(2);
1376  const double arcPosXS = arcPosFunc.calcArcPos(xS);
1377 
1378  FloatArray xE = line.giveVertex(2);
1379  xE.resizeWithValues(2);
1380  const double arcPosXE = arcPosFunc.calcArcPos(xE);
1381 
1382  if ( el->mInteriorSegmentsPointsFine.size() == 0 ) {
1383  Line newLine(xS, xE);
1384  el->mInteriorSegmentsFine.push_back(newLine);
1385  } else {
1386  while ( i < el->mInteriorSegmentsPointsFine.size() ) {
1387 
1388  const FloatArray &x = el->mInteriorSegmentsPointsFine[i];
1389  const double arcPosX = arcPosFunc.calcArcPos(x);
1390 
1391  if ( arcPosX < arcPosXS ) {
1392  OOFEM_ERROR("Error in PrescribedGradientBCWeak :: createTractionMesh.")
1393  }
1394 
1395  if ( arcPosX < arcPosXE ) {
1396  // Split from start pos to x
1397  Line newLine(xS, x);
1398  el->mInteriorSegmentsFine.push_back(newLine);
1399 
1400  xS = x;
1401  } else {
1402  // Split from x to end pos
1403  Line newLine(xS, xE);
1404  el->mInteriorSegmentsFine.push_back(newLine);
1405 
1406  break;
1407  }
1408 
1409  if ( i == (el->mInteriorSegmentsPointsFine.size()-1) ) {
1410  // Split from x to end pos
1411  Line newLine(xS, xE);
1412  el->mInteriorSegmentsFine.push_back(newLine);
1413  }
1414 
1415  i++;
1416  }
1417  }
1418  }
1419  }
1420 
1421  // Create integration rules
1422  for ( TracSegArray* el : mpTracElNew ) {
1423  el->setupIntegrationRuleOnEl();
1424  }
1425 
1426  // Write discontinuity points to debug vtk
1427  std :: vector< FloatArray > discPoints;
1428  for ( TracSegArray* el : mpTracElNew ) {
1429 
1430  discPoints.push_back( el->mInteriorSegments[0].giveVertex(1) );
1431  discPoints.push_back( el->mInteriorSegments.back().giveVertex(2) );
1432  }
1433 
1434 // std :: string fileName("DiscontPoints.vtk");
1435 // XFEMDebugTools :: WritePointsToVTK(fileName, discPoints);
1436 
1437 }
1438 
1439 void PrescribedGradientBCWeak :: splitSegments(std :: vector< TracSegArray * > &ioElArray)
1440 {
1441  std :: vector< TracSegArray * > newArray;
1442 
1443  for ( TracSegArray *el : ioElArray ) {
1444 
1445  for ( Line line : el->mInteriorSegments ) {
1446  TracSegArray *newEl = new TracSegArray();
1447  newEl->mInteriorSegments.push_back(line);
1448  newArray.push_back(newEl);
1449  }
1450 
1451  delete el;
1452  }
1453 
1454  ioElArray = std::move(newArray);
1455 }
1456 
1458 {
1459  TimeStep *tStep = this->giveDomain()->giveEngngModel()->giveCurrentStep();
1460 
1461  double maxDamage = 0.0, damageTol = 0.01;
1462  for ( auto &gp: *el->giveDefaultIntegrationRulePtr() ) {
1463  FloatArray damage;
1464  el->giveIPValue(damage, gp, IST_DamageScalar, tStep);
1465 // printf("damage: "); damage.printYourself();
1466  maxDamage = std::max(maxDamage, damage(0));
1467  }
1468 
1469 // printf("maxDamage: %e\n", maxDamage);
1470 
1471  return maxDamage >= damageTol;
1472 }
1473 
1474 
1475 void PrescribedGradientBCWeak :: assembleTangentGPContributionNew(FloatMatrix &oTangent, TracSegArray &iEl, GaussPoint &iGP, const double &iScaleFactor, const FloatArray &iBndCoord)
1476 {
1478  double detJ = 0.5 * iEl.giveLength();
1479 
1481  // Compute traction N-matrix
1482  // For now, assume piecewise constant approx
1483  FloatArray Ntrac = FloatArray { 1.0*mTracDofScaling };
1484  FloatMatrix NtracMat;
1485  NtracMat.beNMatrixOf(Ntrac, dim);
1486 
1487 
1489  // Compute displacement N-matrix
1490  // Identify the displacement element
1491  // we are currently standing in
1492  // and compute local coordinates on
1493  // the displacement element
1495  FloatArray dispElLocCoord, closestPoint;
1496  Element *dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, iBndCoord );
1497 
1498 
1499  // Compute basis functions
1500  XfemElementInterface *xfemElInt = dynamic_cast< XfemElementInterface * >( dispEl );
1501  FloatMatrix NdispMat;
1502 
1503  if ( xfemElInt != NULL && domain->hasXfemManager() ) {
1504  // If the element is an XFEM element, we use the XfemElementInterface to compute the N-matrix
1505  // of the enriched element.
1506  xfemElInt->XfemElementInterface_createEnrNmatrixAt(NdispMat, dispElLocCoord, * dispEl, false);
1507  } else {
1508  // Otherwise, use the usual N-matrix.
1509  const int numNodes = dispEl->giveNumberOfDofManagers();
1510  FloatArray N(numNodes);
1511 
1512  const int dim = dispEl->giveSpatialDimension();
1513 
1514  NdispMat.resize(dim, dim * numNodes);
1515  NdispMat.zero();
1516  dispEl->giveInterpolation()->evalN( N, dispElLocCoord, FEIElementGeometryWrapper(dispEl) );
1517 
1518  NdispMat.beNMatrixOf(N, dim);
1519  }
1520 
1521  FloatMatrix contrib;
1522  contrib.beTProductOf(NtracMat, NdispMat);
1523  contrib.times( iScaleFactor * detJ * iGP.giveWeight() );
1524 
1525  oTangent = contrib;
1526 }
1527 
1529 {
1530  const double distTol = 1.0e-12;
1531 
1532  if ( iPos [ 0 ] > mUC [ 0 ] - distTol ) {
1533  return true;
1534  }
1535 
1536  if ( iPos [ 1 ] > mUC [ 1 ] - distTol ) {
1537  return true;
1538  }
1539 
1540  return false;
1541 }
1542 
1544 {
1545 //#if 0
1546  if ( mMirrorFunction == 0 ) {
1547  oPosMinus = iPosPlus;
1548  const double distTol = 1.0e-12;
1549 
1550 // if ( iPosPlus.distance(mUC) < distTol ) {
1551 // printf("iPosPlus: %.12e %.12e\n", iPosPlus [ 0 ], iPosPlus [ 1 ]);
1552 // OOFEM_ERROR("Unmappable point.")
1553 // }
1554 
1555  bool mappingPerformed = false;
1556 
1557  if ( iPosPlus [ 0 ] > mUC [ 0 ] - distTol ) {
1558  oPosMinus [ 0 ] = mLC [ 0 ];
1559  mappingPerformed = true;
1560  return;
1561  }
1562 
1563  if ( iPosPlus [ 1 ] > mUC [ 1 ] - distTol ) {
1564  oPosMinus [ 1 ] = mLC [ 1 ];
1565  mappingPerformed = true;
1566  return;
1567  }
1568 
1569  if ( !mappingPerformed ) {
1570  iPosPlus.printYourself();
1571  OOFEM_ERROR("Mapping failed.")
1572  }
1573 
1574  // printf("iPosPlus: "); iPosPlus.printYourself();
1575  // printf("oPosMinus: "); oPosMinus.printYourself();
1576  //#else
1577  }
1578  else {
1579 
1580 #if 1
1581 
1582  const double distTol = 1.0e-12;
1583 // bool mappingPerformed = false;
1584 
1586  FloatArray t = {n(1),-n(0)};
1587  t.normalize();
1588 
1589  double l_s = mUC[0] - mLC[0];
1590 
1591  // Compute angle
1592  double alpha = 0.0, a = 0.0;
1593  if( fabs(t(0)) > 1.0e-6 && fabs(t(1)) > 1.0e-6 ) {
1594  alpha = atan(t(1)/t(0));
1595 
1596  if ( alpha > 45.0*M_PI/180.0 ) {
1597  a = l_s/tan(alpha);
1598  } else {
1599  a = l_s*tan(alpha);
1600  }
1601  } else {
1602  // 90 degrees or 0 degrees
1603 // alpha = 1.57079632679490e+00;
1604  a = 0.0;
1605  }
1606 
1607  if ( alpha > 45.0*M_PI/180.0 ) {
1608 
1609  // alpha > 45 degrees
1610 
1611  if ( iPosPlus [ 0 ] > mUC [ 0 ] - distTol ) {
1612  // Gamma_1_plus
1613  oPosMinus = {0.0, iPosPlus[1]};
1614  return;
1615  }
1616 
1617  if ( iPosPlus [ 1 ] > mUC [ 1 ] - distTol ) {
1618  // Gamma_2_plus
1619 
1620  if ( iPosPlus[0] < a ) {
1621  oPosMinus = {l_s - a + iPosPlus[0], 0.0};
1622  return;
1623  } else {
1624  oPosMinus = {iPosPlus[0] - a, 0.0};
1625  return;
1626  }
1627 
1628  }
1629 
1630  } else {
1631 
1632  // alpha <= 45 degrees
1633 
1634  if ( iPosPlus [ 0 ] > mUC [ 0 ] - distTol ) {
1635  // Gamma_1_plus
1636  if ( iPosPlus[1] < a ) {
1637  oPosMinus = {0.0, l_s - a + iPosPlus[1]};
1638  return;
1639  } else {
1640  oPosMinus = {0.0, iPosPlus[1] - a};
1641  return;
1642  }
1643 
1644  }
1645 
1646  if ( iPosPlus [ 1 ] > mUC [ 1 ] - distTol ) {
1647  // Gamma_2_plus
1648 
1649  oPosMinus = {iPosPlus[0], 0.0};
1650  return;
1651  }
1652 
1653  }
1654 
1655 #else
1656 
1657 #endif
1658  }
1659  //#endif
1660 }
1661 
1663 {
1664 //#if 0
1665  if ( mMirrorFunction == 0 ) {
1666 
1667  const double l_box = mUC(0) - mLC(0);
1668 
1669  oPosPlus = iPosMinus;
1670 // const double distTol = 1.0e-16;
1671  const double distTol = l_box*1.0e-10;
1672 
1673 // if ( iPosMinus.distance(mLC) < distTol ) {
1674 // printf("iPosMinus: %.12e %.12e\n", iPosMinus [ 0 ], iPosMinus [ 1 ]);
1675 // OOFEM_ERROR("Unmappable point.")
1676 // }
1677 
1678  double mappingPerformed = false;
1679 
1680  if ( iPosMinus [ 0 ] < mLC [ 0 ] + distTol ) {
1681  oPosPlus [ 0 ] = mUC [ 0 ];
1682  mappingPerformed = true;
1683  return;
1684  }
1685 
1686  if ( iPosMinus [ 1 ] < mLC [ 1 ] + distTol ) {
1687  oPosPlus [ 1 ] = mUC [ 1 ];
1688  mappingPerformed = true;
1689  return;
1690  }
1691 
1692  if ( !mappingPerformed ) {
1693  iPosMinus.printYourself();
1694  OOFEM_ERROR("Mapping failed.")
1695  }
1696 //#else
1697  } else {
1698 
1699 #if 1
1700 
1701  const double distTol = 1.0e-12;
1702 // bool mappingPerformed = false;
1703 
1705  FloatArray t = {n(1),-n(0)};
1706  t.normalize();
1707 
1708  double l_s = mUC[0] - mLC[0];
1709 
1710  // Compute angle
1711  double alpha = 0.0, a = 0.0;
1712  if ( fabs(t(0)) > 1.0e-6 && fabs(t(1)) > 1.0e-6 ) {
1713  alpha = atan(t(1)/t(0));
1714 
1715  if ( alpha > 45.0*M_PI/180.0 ) {
1716  a = l_s/tan(alpha);
1717  } else {
1718  a = l_s*tan(alpha);
1719  }
1720  } else {
1721  // 90 degrees
1722  a = 0.0;
1723  }
1724 
1725 // printf("t(1)/t(0): %e\n", t(1)/t(0));
1726 // printf("a: %e\n", a);
1727 
1728  if ( alpha > 45.0*M_PI/180.0 ) {
1729 
1730  // alpha > 45 degrees
1731 
1732  if ( iPosMinus [ 0 ] < mLC [ 0 ] + distTol ) {
1733  // Gamma_1_minus
1734  oPosPlus = {l_s, iPosMinus[1]};
1735  return;
1736  }
1737 
1738  if ( iPosMinus [ 1 ] < mLC [ 1 ] + distTol ) {
1739  // Gamma_2_minus
1740 
1741  if ( iPosMinus[0] < l_s - a) {
1742  oPosPlus = {iPosMinus[0] + a, l_s};
1743  return;
1744  }
1745  else {
1746  oPosPlus = {iPosMinus[0] - (l_s - a), l_s};
1747  return;
1748  }
1749 
1750  }
1751 
1752  } else {
1753  // alpha <= 45 degrees
1754 
1755  if ( iPosMinus [ 0 ] < mLC [ 0 ] + distTol ) {
1756  // Gamma_1_minus
1757 
1758  if ( iPosMinus[1] < l_s - a ) {
1759  oPosPlus = {l_s, iPosMinus[1] + a};
1760  return;
1761  } else {
1762  oPosPlus = {l_s, iPosMinus[1] - (l_s - a) };
1763  return;
1764  }
1765  }
1766 
1767  if ( iPosMinus [ 1 ] < mLC [ 1 ] + distTol ) {
1768  // Gamma_2_minus
1769 
1770  oPosPlus = {iPosMinus[0], l_s};
1771  return;
1772 
1773  }
1774 
1775  }
1776 #endif
1777  }
1778 //#endif
1779 }
1780 
1781 
1783 {
1784  // Compute LC and UC by assuming a rectangular domain.
1785  int numNodes = iDomain.giveNumberOfDofManagers();
1786  int nsd = iDomain.giveNumberOfSpatialDimensions();
1787 
1788  FloatArray lc = * ( iDomain.giveDofManager(1)->giveCoordinates() );
1789  FloatArray uc = * ( iDomain.giveDofManager(1)->giveCoordinates() );
1790 
1791  for ( int i = 1; i <= numNodes; i++ ) {
1792  DofManager *dMan = iDomain.giveDofManager(i);
1793  const FloatArray &coord = * ( dMan->giveCoordinates() );
1794 
1795  for ( int j = 0; j < nsd; j++ ) {
1796  if ( coord [ j ] < lc [ j ] ) {
1797  lc [ j ] = coord [ j ];
1798  }
1799 
1800  if ( coord [ j ] > uc [ j ] ) {
1801  uc [ j ] = coord [ j ];
1802  }
1803  }
1804  }
1805 
1806  oLC = std::move(lc);
1807  oUC = std::move(uc);
1808 }
1809 
1811 {
1812  const double distTol = 1.0e-12;
1813 
1814  if ( iPos [ 0 ] > mUC [ 0 ] - distTol ) {
1815  return 0;
1816  }
1817 
1818  if ( iPos [ 1 ] > mUC [ 1 ] - distTol ) {
1819  return 1;
1820  }
1821 
1822  if ( iPos [ 0 ] < mLC [ 0 ] + distTol ) {
1823  return 2;
1824  }
1825 
1826  if ( iPos [ 1 ] < mLC [ 1 ] + distTol ) {
1827  return 3;
1828  }
1829 
1830  OOFEM_ERROR("Could not identify side index.")
1831 
1832  return -1;
1833 }
1834 
1835 void PrescribedGradientBCWeak :: findHoleCoord(std::vector<FloatArray> &oHoleCoordUnsorted, std::vector<FloatArray> &oAllCoordUnsorted)
1836 {
1837  Set *setPointer = this->giveDomain()->giveSet(this->set);
1838  const IntArray &boundaries = setPointer->giveBoundaryList();
1839  IntArray bNodes;
1840 
1841  // Loop over boundary nodes and check how many times they occur:
1842  // 1 -> at the edge of an inclusion, therefore must be retained
1843  // 2 -> connected to two segments, optional to keep
1844 
1845  std::unordered_map<int,int> map_bnd_node_ind_to_num_occurences;
1846  for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
1847 
1848  int elIndex = boundaries.at(pos * 2 - 1);
1849  Element *e = this->giveDomain()->giveElement( elIndex );
1850  int boundary = boundaries.at(pos * 2);
1851  e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
1852  DofManager *startNode = e->giveDofManager(bNodes [ 0 ]);
1853  int startNodeInd = startNode->giveNumber();
1854  DofManager *endNode = e->giveDofManager(bNodes [ 1 ]);
1855  int endNodeInd = endNode->giveNumber();
1856 
1857 
1858  auto res = map_bnd_node_ind_to_num_occurences.find(startNodeInd);
1859  if ( res != map_bnd_node_ind_to_num_occurences.end() ) {
1860  map_bnd_node_ind_to_num_occurences[startNodeInd]++;
1861  } else {
1862  map_bnd_node_ind_to_num_occurences[startNodeInd] = 1;
1863  }
1864 
1865  res = map_bnd_node_ind_to_num_occurences.find(endNodeInd);
1866  if ( res != map_bnd_node_ind_to_num_occurences.end() ) {
1867  map_bnd_node_ind_to_num_occurences[endNodeInd]++;
1868  } else {
1869  map_bnd_node_ind_to_num_occurences[endNodeInd] = 1;
1870  }
1871 
1872  }
1873 
1874 
1875  for ( auto it = map_bnd_node_ind_to_num_occurences.begin(); it != map_bnd_node_ind_to_num_occurences.end(); ++it ) {
1876 
1877  bool mandatory_to_keep = false;
1878  if ( it->second == 1 ) {
1879  mandatory_to_keep = true;
1880  }
1881 
1882  DofManager *bndNode = domain->giveDofManager(it->first);
1883  const FloatArray &x = * ( bndNode->giveCoordinates() );
1884  FloatArray xPlus = x;
1885 
1886  if ( !boundaryPointIsOnActiveBoundary(x) ) {
1887  giveMirroredPointOnGammaPlus(xPlus, x);
1888  }
1889 
1890  if ( mandatory_to_keep ) {
1891  oHoleCoordUnsorted.push_back(xPlus);
1892  }
1893 
1894  oAllCoordUnsorted.push_back(xPlus);
1895  }
1896 
1897 
1898 }
1899 
1900 void PrescribedGradientBCWeak :: findCrackBndIntersecCoord(std::vector<FloatArray> &oHoleCoordUnsorted)
1901 {
1902  Set *setPointer = this->giveDomain()->giveSet(this->set);
1903  const IntArray &boundaries = setPointer->giveBoundaryList();
1904  IntArray bNodes;
1905 
1906  for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
1907  Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
1908  int boundary = boundaries.at(pos * 2);
1909 
1910  e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
1911 
1912  // Add the start and end nodes of the segment
1913  DofManager *startNode = e->giveDofManager(bNodes [ 0 ]);
1914  const FloatArray &xS = * ( startNode->giveCoordinates() );
1915 
1916  DofManager *endNode = e->giveDofManager(bNodes [ 1 ]);
1917  const FloatArray &xE = * ( endNode->giveCoordinates() );
1918 
1919  FloatArray xC;
1920  xC.beScaled(0.5, xS);
1921  xC.add(0.5, xE);
1922 
1923 
1924  // Add the points where cracks intersect the boundary
1925  XfemElementInterface *xfemElInt = dynamic_cast< XfemElementInterface * >( e );
1926 
1927  if ( xfemElInt != NULL && domain->hasXfemManager() ) {
1928  std :: vector< Line >segments;
1929  std :: vector< FloatArray >intersecPoints;
1930  xfemElInt->partitionEdgeSegment(boundary, segments, intersecPoints, mTangDistPadding);
1931 
1932  for ( auto x : intersecPoints ) {
1933 
1934  if ( pointIsOnGammaPlus(x) ) {
1935  oHoleCoordUnsorted.push_back(x);
1936  } else {
1937  FloatArray xPlus;
1938  giveMirroredPointOnGammaPlus(xPlus, x);
1939  oHoleCoordUnsorted.push_back( std::move(xPlus) );
1940  }
1941  }
1942 
1943  }
1944 
1945  }
1946 
1947 
1948  // Also add traction nodes where cohesive zone elements intersect the boundary
1949  // TODO: Check later, this should already be covered by the new approach for
1950  // identifying inclusions.
1951 }
1952 
1953 void PrescribedGradientBCWeak :: findPeriodicityCoord(std::vector<FloatArray> &oHoleCoordUnsorted)
1954 {
1955  const double l_s = mUC[0] - mLC[0];
1956 
1958  FloatArray t = {n(1),-n(0)};
1959 
1960  if ( mMirrorFunction == 1 || mMirrorFunction == 2 ) {
1961 
1962  if ( fabs(n(1)) <= fabs(n(0)) ) {
1963  // a <= l_s/2
1964  double a = 0.5*l_s*( 1.0 + n(1)/n(0) );
1965 
1966  FloatArray p1 = {2.0*a, 0.0};
1967  FloatArray p1Plus;
1968  giveMirroredPointOnGammaPlus(p1Plus, p1);
1969  oHoleCoordUnsorted.push_back(std::move(p1Plus));
1970  } else {
1971  // a > l_s/2
1972  double c = l_s - 0.5*l_s*( 1.0 + t(1)/t(0) );
1973 
1974  FloatArray p1 = {l_s, l_s-2.0*c};
1975  oHoleCoordUnsorted.push_back(std::move(p1));
1976 
1977  FloatArray p2 = {0, 2.0*c};
1978  FloatArray p2Plus;
1979  giveMirroredPointOnGammaPlus(p2Plus, p2);
1980  oHoleCoordUnsorted.push_back(std::move(p2Plus));
1981  }
1982  }
1983 
1984 }
1985 
1986 void PrescribedGradientBCWeak :: removeClosePoints(std::vector<FloatArray> &ioCoords, const double &iAbsTol)
1987 {
1988  if ( ioCoords.size() == 0 ) {
1989  return;
1990  }
1991 
1992  const double tol2 = iAbsTol*iAbsTol;
1993 
1994  std::vector<FloatArray> tmp;
1995  tmp.push_back(ioCoords[0]);
1996  size_t j = 0;
1997 
1998  for ( size_t i = 1; i < ioCoords.size(); i++ ) {
1999  if ( ioCoords[i].distance_square(tmp[j]) > tol2 ) {
2000  tmp.push_back(ioCoords[i]);
2001  j++;
2002  }
2003  }
2004 
2005  ioCoords = std::move(tmp);
2006 }
2007 
2009 {
2010  // Idea: Loop over segments and check if the center point of each
2011  // segment is located on an element. If not, the segment is
2012  // located over a hole and needs to be removed.
2013 
2014  std :: vector< Line > tmp;
2015 
2017  const double tol2 = iAbsTol*iAbsTol;
2018 
2019  for ( auto &l : ioTSeg.mInteriorSegments ) {
2020  const FloatArray &xS = l.giveVertex(1);
2021  const FloatArray &xE = l.giveVertex(2);
2022  FloatArray xPlus = {0.5*(xS[0]+xE[0]), 0.5*(xS[1]+xE[1])};
2023 
2024  FloatArray lcoordsPlus, closestPlus;
2025  localizer->giveElementClosestToPoint(lcoordsPlus, closestPlus, xPlus);
2026 
2027  FloatArray xMinus;
2028  giveMirroredPointOnGammaMinus(xMinus, xPlus);
2029  FloatArray lcoordsMinus, closestMinus;
2030  localizer->giveElementClosestToPoint(lcoordsMinus, closestMinus, xMinus);
2031 
2032  if ( !(xPlus.distance_square(closestPlus) > tol2 || xMinus.distance_square(closestMinus) > tol2) ) {
2033  tmp.push_back(l);
2034  }
2035  }
2036 
2037  ioTSeg.mInteriorSegments = std::move(tmp);
2038 }
2039 
2040 } /* namespace oofem */
bool contains(int value) const
Definition: intarray.h:283
The base class for all spatial localizers.
void findCrackBndIntersecCoord(std::vector< FloatArray > &oHoleCoordUnsorted)
void setField(int item, InputFieldType id)
bool mMeshIsPeriodic
true -> the traction lives only on gammaPlus, so that we get strong periodicity as a special case...
The representation of EngngModel default unknown numbering.
const FloatArray & giveVertex(int n) const
Definition: geometry.h:124
std::vector< Line > mInteriorSegments
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in full form.
Definition: element.C:1257
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)=0
Evaluates the array of interpolation functions (shape functions) at given point.
Class and object Domain.
Definition: domain.h:115
virtual int assemble(const IntArray &loc, const FloatMatrix &mat)=0
Assembles sparse matrix from contribution of local elements.
virtual IntegrationRule * giveDefaultIntegrationRulePtr()
Access method for default integration rule.
Definition: element.h:822
virtual int giveNumberOfInternalDofManagers()
Gives the number of internal dof managers.
void findHoleCoord(std::vector< FloatArray > &oHoleCoordUnsorted, std::vector< FloatArray > &oAllCoordUnsorted)
virtual void assembleVector(FloatArray &answer, TimeStep *tStep, CharType type, ValueModeType mode, const UnknownNumberingScheme &s, FloatArray *eNorm=NULL)
Assembles B.C.
virtual void giveInputRecord(DynamicInputRecord &input)
void giveTractionElCoord(size_t iElInd, FloatArray &oStartCoord, FloatArray &oEndCoord) const
int giveGlobalNumber() const
Definition: dofmanager.h:501
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition: femcmpnn.h:82
int giveNumberOfDofManagers() const
Returns number of dof managers in domain.
Definition: domain.h:432
Base class for all matrices stored in sparse format.
Definition: sparsemtrx.h:60
IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: activebc.h:75
SparseMtrx * createSparseMtrx(SparseMtrxType type)
Creates new instance of sparse matrix corresponding to given keyword.
Definition: classfactory.C:105
const IntArray & giveBoundaryList()
Returns list of element boundaries within set.
Definition: set.C:140
FloatArray mUC
Upper corner of domain (assuming a rectangular RVE)
Class for homogenization of applied gradients.
Provides Xfem interface for an element.
std::unique_ptr< IntegrationRule > mIntRule
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
ValueModeType
Type representing the mode of UnknownType or CharType, or similar types.
Definition: valuemodetype.h:78
int mTractionNodeSpacing
Use every (mTractionNodeSpacing) displacement nodes when constructing the traction element mesh...
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
#define _IFT_PrescribedGradientBCWeak_TractionInterpOrder
Solves a static structural problem.
virtual FloatArray * giveCoordinates()
Definition: dofmanager.h:382
EngngModel * giveEngngModel()
Returns engineering model to which receiver is associated.
Definition: domain.C:433
void findPeriodicityCoord(std::vector< FloatArray > &oHoleCoordUnsorted)
Implementation for assembling tangent matrices in standard monolithic FE-problems.
int giveSideIndex(const FloatArray &iPos) const
void beVectorForm(const FloatMatrix &aMatrix)
Reciever will be a vector with 9 components formed from a 3x3 matrix.
Definition: floatarray.C:992
DiscontinuousSegmentIntegrationRule provides integration over a discontinuous boundary segment...
virtual void assembleExtraDisplock(SparseMtrx &answer, TimeStep *tStep, CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)
double giveTargetTime()
Returns target time.
Definition: timestep.h:146
Abstract base class for all finite elements.
Definition: element.h:145
void giveTraction(size_t iElInd, FloatArray &oStartTraction, FloatArray &oEndTraction, ValueModeType mode, TimeStep *tStep)
Node * mpDisplacementLock
Lock displacements in one node if periodic.
Base class for dof managers.
Definition: dofmanager.h:113
#define S(p)
Definition: mdm.C:481
int giveNumber()
Returns domain number.
Definition: domain.h:266
int giveNumberOfSpatialDimensions()
Returns number of spatial dimensions.
Definition: domain.C:1067
virtual void giveInputRecord(DynamicInputRecord &input)
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
int mTractionInterpOrder
Order of interpolation for traction (0->piecewise constant, 1->piecewise linear)
bool pointIsOnGammaPlus(const FloatArray &iPos) const
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
XfemManager * giveXfemManager()
Definition: domain.C:375
void sort(IntArray &arry, operation op)
Sorts the receiver using quicksort algorithm.
Definition: intarray.h:416
virtual int giveNumberOfDofManagers() const
Definition: element.h:656
virtual FEInterpolation * giveInterpolation() const
Definition: element.h:629
#define _IFT_PrescribedGradientBCWeak_NumTractionNodeSpacing
void computeIntForceGPContrib(FloatArray &oContrib_disp, IntArray &oDisp_loc_array, FloatArray &oContrib_trac, IntArray &oTrac_loc_array, TracSegArray &iEl, GaussPoint &iGP, int iDim, TimeStep *tStep, const FloatArray &iBndCoord, const double &iScaleFac, ValueModeType mode, CharType type, const UnknownNumberingScheme &s)
int mMirrorFunction
Mirror function (i.e.
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
void giveTractionElNormal(size_t iElInd, FloatArray &oNormal, FloatArray &oTangent) const
virtual void computeField(FloatArray &sigma, TimeStep *tStep)
Computes the homogenized, macroscopic, field (stress).
virtual void assemble(SparseMtrx &answer, TimeStep *tStep, const MatrixAssembler &ma, const UnknownNumberingScheme &s, Domain *domain)
Assembles characteristic matrix of required type into given sparse matrix.
Definition: engngm.C:803
Class representing "master" degree of freedom.
Definition: masterdof.h:92
void giveTractionLocationArray(IntArray &rows, CharType type, const UnknownNumberingScheme &s)
void beScaled(double s, const FloatArray &b)
Sets receiver to be .
Definition: floatarray.C:146
#define max
void giveMirroredPointOnGammaMinus(FloatArray &oPosMinus, const FloatArray &iPosPlus) const
std::unique_ptr< SparseMtrx > stiffnessMatrix
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
void giveUnknownVector(FloatArray &answer, const IntArray &dofMask, ValueModeType mode, TimeStep *tStep, bool padding=false)
Assembles the vector of unknowns in global c.s for given dofs of receiver.
Definition: dofmanager.C:685
void giveBoundaries(IntArray &oBoundaries)
const IntArray & giveTracDofIDs() const
#define E(p)
Definition: mdm.C:368
Element * giveElement(int n)
Service for accessing particular domain fe element.
Definition: domain.C:160
#define M_PI
Definition: mathfem.h:52
#define OOFEM_ERROR(...)
Definition: error.h:61
void removeClosePoints(std::vector< FloatArray > &ioCoords, const double &iAbsTol)
void createTractionMesh(bool iEnforceCornerPeriodicity, int iNumSides)
void clear()
Clears the array (zero size).
Definition: intarray.h:177
DofIDItem
Type representing particular dof type.
Definition: dofiditem.h:86
const IntArray & giveDispLockDofIDs() const
virtual void giveTractionLocationArray(IntArray &rows, const UnknownNumberingScheme &s)
SparseMtrxType
Enumerative type used to identify the sparse matrix type.
virtual Node * giveNodeClosestToPoint(const FloatArray &coords, double maxDist)=0
Returns the node closest to the given coordinate.
virtual double giveWeight()
Returns integration weight of receiver.
Definition: gausspoint.h:181
double calcArcPos(const FloatArray &iPos) const
void resizeWithData(int, int)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1369
void giveLocationArray(IntArray &locationArray, const UnknownNumberingScheme &s, IntArray *dofIds=NULL) const
Returns the location array (array of code numbers) of receiver for given numbering scheme...
Definition: element.C:390
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
Set of elements, boundaries, edges and/or nodes.
Definition: set.h:66
SpatialLocalizer * giveSpatialLocalizer()
Returns receiver&#39;s associated spatial localizer.
Definition: domain.C:1184
double mTangDistPadding
Parameter for creation of traction mesh.
Abstract base class allowing to control the way, how equations are assigned to individual DOFs...
void computeExtForceElContrib(FloatArray &oContrib, TracSegArray &iEl, int iDim, TimeStep *tStep)
#define _IFT_PrescribedGradientBCWeak_TangDistPadding
void XfemElementInterface_createEnrNmatrixAt(FloatMatrix &oAnswer, const FloatArray &iLocCoord, Element &iEl, bool iSetDiscontContribToZero)
Creates enriched N-matrix.
virtual void giveLocationArrays(std::vector< IntArray > &rows, std::vector< IntArray > &cols, CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)
Gives a list of location arrays that will be assembled.
Wrapper around element definition to provide FEICellGeometry interface.
Definition: feinterpol.h:95
Set * giveSet(int n)
Service for accessing particular domain set.
Definition: domain.C:363
FloatArray mPeriodicityNormal
Periodicity direction.
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
#define N(p, q)
Definition: mdm.C:367
#define _IFT_PrescribedGradientBCWeak_DuplicateCornerNodes
Abstract base class for all active boundary conditions.
Definition: activebc.h:63
double calcArcPos(const FloatArray &iPos) const
static int giveVoigtVectorMask(IntArray &answer, MaterialMode mmode)
Returns a mask of the vector indicies corresponding to components in a general (non-symmetric) second...
void compute_x_times_N_1(FloatMatrix &o_x_times_N)
void resizeWithValues(int s, int allocChunk=0)
Checks size of receiver towards requested bounds.
Definition: floatarray.C:615
std::vector< Line > mInteriorSegmentsFine
void beTProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix^T and anArray.
Definition: floatarray.C:708
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
double distance_square(const FloatArray &iP1, const FloatArray &iP2, double &oXi, double &oXiUnbounded) const
Definition: floatarray.C:499
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
int mNumTractionNodesAtIntersections
If traction nodes should be inserted where cracks intersect the RVE boundary.
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
void partitionEdgeSegment(int iBndIndex, std::vector< Line > &oSegments, std::vector< FloatArray > &oIntersectionPoints, const double &iTangDistPadding=0.0)
Partition a boundary segment to account for cracks cutting the boundary.
virtual void scale(double s)
Scales the receiver according to given value.
int giveDofManPlaceInArray(int iGlobalDofManNum) const
Returns the array index of the dofman with global number iGlobalDofManNum, so that it can be fetched ...
Definition: domain.C:196
const FloatArray & giveGlobalCoordinates()
Definition: gausspoint.h:160
void appendDof(Dof *dof)
Adds the given Dof into the receiver.
Definition: dofmanager.C:134
void beNMatrixOf(const FloatArray &n, int nsd)
Assigns the receiver to be a repeated diagonal matrix.
Definition: floatmatrix.C:505
void giveMirroredPointOnGammaPlus(FloatArray &oPosPlus, const FloatArray &iPosMinus) const
bool hasXfemManager()
Definition: domain.C:386
#define _IFT_PrescribedGradientBCWeak_TracDofScaling
virtual void boundaryGiveNodes(IntArray &answer, int boundary)=0
Gives the boundary nodes for requested boundary number.
Class representing vector of real numbers.
Definition: floatarray.h:82
SparseLinearSystemNM * createSparseLinSolver(LinSystSolverType st, Domain *d, EngngModel *m)
Creates new instance of SparseLinearSystemNM corresponding to given type.
Definition: classfactory.C:120
#define _IFT_PrescribedGradientBCWeak_MirrorFunction
void assembleTangentGPContributionNew(FloatMatrix &oTangent, TracSegArray &iEl, GaussPoint &iGP, const double &iScaleFactor, const FloatArray &iBndCoord)
virtual int giveSpatialDimension()
Returns the element spatial dimension (1, 2, or 3).
Definition: element.C:1347
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
This class manages the xfem part.
Definition: xfemmanager.h:109
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
void assemble(const FloatArray &fe, const IntArray &loc)
Assembles the array fe (typically, the load vector of a finite element) into the receiver, using loc as location array.
Definition: floatarray.C:551
virtual DofManager * giveInternalDofManager(int i)
Gives an internal dof manager from receiver.
virtual void giveDisplacementLocationArray(IntArray &rows, const UnknownNumberingScheme &s)
CharType
Definition: chartype.h:87
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
virtual Element * giveElementClosestToPoint(FloatArray &lcoords, FloatArray &closest, const FloatArray &coords, int region=0)=0
Returns the element closest to a given point.
Class representing the general Input Record.
Definition: inputrecord.h:101
virtual void printYourself() const
Print receiver on stdout.
Definition: floatarray.C:747
std::unique_ptr< Node > mFirstNode
virtual void postInitialize()
Performs post initialization steps.
void append(const FloatArray &a)
Appends array to reciever.
Definition: floatarray.C:664
FloatArray mLC
Lower corner of domain (assuming a rectangular RVE)
virtual void assemble(SparseMtrx &answer, TimeStep *tStep, CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s, double scale=1.0)
Assembles B.C.
void followedBy(const IntArray &b, int allocChunk=0)
Appends array b at the end of receiver.
Definition: intarray.C:145
void computeDomainBoundingBox(Domain &iDomain, FloatArray &oLC, FloatArray &oUC)
Class implementing single timer, providing wall clock and user time capabilities. ...
Definition: timer.h:46
void compute_x_times_N_2(FloatMatrix &o_x_times_N)
void beDyadicProductOf(const FloatArray &vec1, const FloatArray &vec2)
Assigns to the receiver the dyadic product .
Definition: floatmatrix.C:492
void beTProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:367
void giveTractionElArcPos(size_t iElInd, double &oXiStart, double &oXiEnd) const
void removeSegOverHoles(TracSegArray &ioTSeg, const double &iAbsTol)
Class representing the a dynamic Input Record.
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
void beTranspositionOf(const FloatMatrix &src)
Assigns to the receiver the transposition of parameter.
Definition: floatmatrix.C:323
ClassFactory & classFactory
Definition: classfactory.C:59
virtual void computeTangent(FloatMatrix &E, TimeStep *tStep)
Computes the macroscopic tangent for homogenization problems through sensitivity analysis.
#define _IFT_PrescribedGradientBCWeak_NumTractionNodesAtIntersections
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
Domain * giveDomain() const
Definition: femcmpnn.h:100
virtual void assembleGPContrib(SparseMtrx &answer, TimeStep *tStep, CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s, TracSegArray &iEl, GaussPoint &iGP, double k)
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
void push_back(const double &iVal)
Add one element.
Definition: floatarray.h:118
void beUnitMatrix()
Sets receiver to unity matrix.
Definition: floatmatrix.C:1332
const IntArray & giveRegularDispDofIDs() const
#define _IFT_PrescribedGradientBCWeak_PeriodicityNormal
Abstract base class representing the "problem" under consideration.
Definition: engngm.h:181
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
virtual TimeStep * giveCurrentStep(bool force=false)
Returns current time step.
Definition: engngm.h:683
int giveSize() const
Definition: intarray.h:203
virtual bool boundaryPointIsOnActiveBoundary(const FloatArray &iPos) const =0
void splitSegments(std::vector< TracSegArray * > &ioElArray)
virtual void giveBoundaryCoordVector(FloatArray &oX, const FloatArray &iPos) const =0
the oofem namespace is to define a context or scope in which all oofem names are defined.
DofManager * giveDofManager(int i) const
Definition: element.C:514
void assemble(const FloatMatrix &src, const IntArray &loc)
Assembles the contribution using localization array into receiver.
Definition: floatmatrix.C:251
const IntArray & giveEnrichedDofIDs() const
Definition: xfemmanager.h:191
void clear()
Sets size of receiver to be an empty matrix. It will have zero rows and zero columns size...
Definition: floatmatrix.h:516
Class implementing node in finite element mesh.
Definition: node.h:87
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
int giveNumber() const
Definition: femcmpnn.h:107
DofManager * giveDofManager(int n)
Service for accessing particular domain dof manager.
Definition: domain.C:314
void negated()
Switches the sign of every coefficient of receiver.
Definition: floatarray.C:739
double normalize()
Normalizes receiver.
Definition: floatarray.C:828
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
std::vector< TracSegArray * > mpTracElNew
Elements for the independent traction discretization.
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
void giveLocationArray(const IntArray &dofIDArry, IntArray &locationArray, const UnknownNumberingScheme &s) const
Returns location array (array containing for each requested dof related equation number) for given nu...
Definition: dofmanager.C:177
void startTimer()
Definition: timer.C:69
virtual double evaluateAtTime(double t)
Returns the value of the function at given time.
Definition: function.C:76
double getUtime()
Returns total user time elapsed in seconds.
Definition: timer.C:105
Class representing integration point in finite element program.
Definition: gausspoint.h:93
Class representing solution step.
Definition: timestep.h:80
void stopTimer()
Definition: timer.C:77
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
bool mDuplicateCornerNodes
0 -> Do not duplicate corner traction nodes 1 -> Duplicate corner traction nodes

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:30 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011