OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
quasicontinuum.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; eitherc
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 
35 #include "../sm/Quasicontinuum/quasicontinuum.h"
36 #include "../sm/EngineeringModels/qclinearstatic.h"
37 #include "qcnode.h"
38 #include "element.h"
39 #include "mathfem.h"
40 
41 #include "classfactory.h"
42 #include "dynamicinputrecord.h"
43 #include "../sm/CrossSections/simplecrosssection.h"
44 #include "../sm/Materials/isolinearelasticmaterial.h"
45 #include "../sm/Materials/anisolinearelasticmaterial.h"
46 #include "interfacetype.h"
47 #include "../sm/Materials/qcmaterialextensioninterface.h"
48 
49 
50 namespace oofem {
51 //REGISTER_Quasicontinuum(QCFullsolveddomain);
52 
54 // Constructor.
55 {}
56 
58 // Destructor
59 {}
60 
61 
62 void
64 // rewrite existing interpolation mesh
65 {
66  int dim = d->giveNumberOfSpatialDimensions();
67  if ( dim == 2 || dim == 3 ) {
68  this->nDimensions = dim;
69  } else {
70  OOFEM_ERROR("Invalid number of dimensions. Only 2d and 3d domains are supported in QC simulation. \n");
71  }
72 }
73 
74 
75 void
76 Quasicontinuum :: setupInterpolationMesh(Domain *d, int generateInterpolationElements, int interpolationElementsMaterialNumber, std :: vector< IntArray > *newMeshNodes)
77 // rewrite existing interpolation mesh
78 {
79  int noNewEl = 0;
80  if ( generateInterpolationElements == 0 ) {
81  // nodes are not needed
82  this->interpolationMeshNodes.clear();
83  // elements are defined in imput file by material number
84  this->interpolationElementNumbers = d->giveElementsWithMaterialNum(interpolationElementsMaterialNumber);
85  } else if ( ( generateInterpolationElements == 1 ) || ( generateInterpolationElements == 2 ) ) {
86  // nodes were (generated and ) loaded from t3d output file
87  this->interpolationMeshNodes = * newMeshNodes;
88  // elements will be placed in the end of element list
90  int noEl = d->giveNumberOfElements();
91  noNewEl = newMeshNodes->size();
92  this->interpolationElementNumbers.resize(noNewEl);
93  for ( int i = 1; i <= noNewEl; i++ ) {
94  this->interpolationElementNumbers.at(i) = noEl + i;
95  }
96  } else {
97  OOFEM_ERROR("Invalid input field \"GenInterpElem\". Value: %d is not supported", generateInterpolationElements);
98  }
99 
101  for ( int i = 1; i <= interpolationElementNumbers.giveSize(); i++ ) {
103  }
104 }
105 
106 
107 
108 void
110 // create new elements. Material will be defined after post initialization
111 {
112  /*
113  * // new crosssection
114  * int ncrosssect = d->giveNumberOfCrossSectionModels() + 1;
115  * d->resizeCrossSectionModels(ncrosssect);
116  * CrossSection *crossSection;
117  * crossSection = classFactory.createCrossSection("SimpleCS", ncrosssect, d);
118  * DynamicInputRecord irCS;
119  * irCS.setField(1.0, _IFT_SimpleCrossSection_thick);
120  * irCS.setField(1.0, _IFT_SimpleCrossSection_width);
121  * crossSection->initializeFrom(&irCS);
122  * d->setCrossSection(ncrosssect, crossSection);
123  */
124 
125  int ninterpelem = interpolationMeshNodes.size();
126 
127  /*
128  * #if DEBUG
129  * if (ninterpelem==0) {
130  * OOFEM_WARNING("No interpolation elements generated");
131  * }
132  * #endif
133  */
134 
135  // add interpolation elements (with new CS and no mat)
136  const char *elemType;
137  if ( interpolationMeshNodes.size() != 0 ) {
138  // for 2d
139  if ( nDimensions == 2 ) {
140  //elemType = "TrPlaneStrain";
141  elemType = "TrPlaneStress2d";
142  //for 3d
143  } else {
144  elemType = "LTRSpace";
145  }
146  }
147 
148  int nelem = d->giveNumberOfElements();
149  int elemNumber;
150  Element *elem;
151  DynamicInputRecord irEl;
152  // over interpolation elements
153  for ( int i = 1; i <= ninterpelem; i++ ) {
154  elemNumber = nelem + i;
155  d->resizeElements(elemNumber);
156  elem = classFactory.createElement(elemType, elemNumber, d);
158  //irEl.setField( ncrosssect, _IFT_Element_crosssect);
159  //irEl.setField( nmat, _IFT_Element_mat);
160  elem->initializeFrom(& irEl);
161  elem->setGlobalNumber(elemNumber);
162  d->setElement(elemNumber, elem);
163  }
164 }
165 
166 void
168 {
169  // new crosssection
170  int ncrosssect = d->giveNumberOfCrossSectionModels() + 1;
171  d->resizeCrossSectionModels(ncrosssect);
172  CrossSection *crossSection;
173  crossSection = classFactory.createCrossSection("SimpleCS", ncrosssect, d);
174  DynamicInputRecord irCS;
177  crossSection->initializeFrom(& irCS);
178  d->setCrossSection(ncrosssect, crossSection);
179 
180  int elNum;
181  for ( int i = 1; i <= interpolationElementNumbers.giveSize(); i++ ) {
182  elNum = interpolationElementNumbers.at(i);
183  d->giveElement(elNum)->setCrossSection(ncrosssect);
184  }
185 }
186 
187 
188 
189 void
191 // Approach 1 --- Hangingnodes (interpolation elementss with zero stiffnes)
192 {
193  // new material with zero stiffness
194  int nmat = d->giveNumberOfMaterialModels() + 1;
195  d->resizeMaterials(nmat);
196  Material *mat;
197  mat = classFactory.createMaterial("IsoLE", nmat, d);
198  DynamicInputRecord irMat;
202  irMat.setField(0.0, _IFT_Material_density);
203  mat->initializeFrom(& irMat);
204  d->setMaterial(nmat, mat);
205 
206 
207  // add material and CS to interpolation elements
208  int ninterpelem = interpolationElementNumbers.giveSize();
209  int elNum;
210  int matNum = nmat;
211  // over interpolation elements
212  for ( int i = 1; i <= ninterpelem; i++ ) {
213  elNum = interpolationElementNumbers.at(i);
214  d->giveElement(elNum)->setMaterial(matNum);
215  }
216 }
217 
218 void
219 Quasicontinuum :: applyApproach2(Domain *d, int homMtrxType, double volumeOfInterpolationMesh)
220 // Approach 2 --- Global homogenization
221 {
222  // decision which links will solved explicitly/homogenized
224  elemList.zero();
226  nodeList.zero();
227 
228  // ovel all elements // TO DO improve by using list of links only
229  for ( int i = 1; i <= d->giveNumberOfElements(); i++ ) {
230  if ( d->giveElement(i)->giveNumberOfNodes() > 2 ) {
231  elemList.at(i) = 1; // all interpolation elements will be active
232  continue; // skip interpolation elements
233  }
234 
235  Element *e = d->giveElement(i);
236  int n1 = e->giveNode(1)->giveQcNodeType();
237  int n2 = e->giveNode(2)->giveQcNodeType();
238  // repnode--repnode - save link to elemList
239  if ( ( n1 == 1 ) && ( n2 == 1 ) ) {
240  elemList.at(i) = 1;
241  nodeList.at( e->giveNode(1)->giveNumber() ) = 1;
242  nodeList.at( e->giveNode(2)->giveNumber() ) = 1;
243 
244  // hanging--hanging
245  } else if ( ( n1 == 2 ) && ( n2 == 2 ) ) {
246  // nothing to store, link will be homogenized-removed
247  // eventualy, link with end/ends outside intero element can be solved explicitly
248 
249  // repnode--hanging
250  } else if ( ( n1 == 1 ) && ( n2 == 2 ) ) {
251  int mn = e->giveNode(1)->giveNumber();
252  int hn = e->giveNode(2)->giveNumber();
253  qcNode *qn = static_cast< qcNode * >( e->giveNode(2) );
254  int masterElem = qn->giveMasterElementNumber();
255  IntArray nodesOfMasterElem = d->giveElement(masterElem)->giveDofManArray();
256  // is masterNode one of nodes of master element?
257  if ( !nodesOfMasterElem.contains(mn) ) { // master node is not vertex of hn master elements
258  elemList.at(i) = 1;
259  nodeList.at(mn) = 1;
260  nodeList.at(hn) = 1;
261  } else {
262  // nothing to store, link will be homogenized-removed
263  }
264 
265  // hanging--repnode
266  } else if ( ( n1 == 2 ) && ( n2 == 1 ) ) {
267  int mn = e->giveNode(2)->giveNumber();
268  int hn = e->giveNode(1)->giveNumber();
269  qcNode *n = static_cast< qcNode * >( e->giveNode(1) );
270  int masterElem = n->giveMasterElementNumber();
271  IntArray nodesOfMasterElem = d->giveElement(masterElem)->giveDofManArray();
272  // is masterNode one of nodes of master element?
273  if ( !nodesOfMasterElem.contains(mn) ) { // master node is not vertex of hn master elements
274  elemList.at(i) = 1;
275  nodeList.at(mn) = 1;
276  nodeList.at(hn) = 1;
277  } else {
278  // nothing to store, link will be homogenized-removed
279  }
280 
281 
282  // no QcNode
283  } else {
284  OOFEM_WARNING("Element %d with non-qcNode can not be homogenized", i);
285  elemList.at(i) = 1;
286  nodeList.at( e->giveNode(1)->giveNumber() ) = 1;
287  nodeList.at( e->giveNode(2)->giveNumber() ) = 1;
288  }
289  } // ovel all links
290 
291 
292  /* // qcRenumbering ir realized instead of this...
293  * // deactivate homogenized nodes
294  * for (int i=1; i<=d->giveNumberOfDofManagers(); i++) {
295  * if ( !nodeList(i) ) {
296  * /// d->giveNode(i)->deactivateYourself();
297  * }
298  * }
299  *
300  * // deactivate homogenized elements
301  * for (int i=1; i<=d->giveNumberOfElements(); i++) {
302  * if ( !elemList(i) ) {
303  * /// d->giveElement(i)->deactivateYourself();
304  * }
305  * }
306  */
307 
308  // set up lists of active nodes and elements
309  QClinearStatic *em = dynamic_cast< QClinearStatic * >( d->giveEngngModel() );
310  if ( em ) {
313  } else {
314  OOFEM_ERROR("Quasicontinuum can be applied only in QClinearStatic engngm");
315  }
316 
317 
318 
319  // global homogenization
320  double S0;
321  FloatMatrix DisoMatrix(6, 6);
322  FloatMatrix *Diso = & DisoMatrix;
323  createGlobalStiffnesMatrix(Diso, S0, d, homMtrxType, volumeOfInterpolationMesh);
324 
325  double homogenizedE, homogenizedNu;
326  homogenizationOfStiffMatrix(homogenizedE, homogenizedNu, Diso);
327 
328 
329  // new material with homogenized stiffness
330  int nmat = d->giveNumberOfMaterialModels() + 1;
331  d->resizeMaterials(nmat);
332  Material *mat;
333  DynamicInputRecord irMat;
334  // isotropic
335  if ( homMtrxType == 1 ) {
336  mat = classFactory.createMaterial("IsoLE", nmat, d);
337  irMat.setField(homogenizedE, _IFT_IsotropicLinearElasticMaterial_e);
338  irMat.setField(homogenizedNu, _IFT_IsotropicLinearElasticMaterial_n);
340  irMat.setField(0.0, _IFT_Material_density);
341  // anisotropic
342  } else if ( homMtrxType == 2 ) {
343  //OOFEM_ERROR("anisotropic homog. is not inmplemented yet");
344  mat = classFactory.createMaterial("AnisoLE", nmat, d);
345  FloatArray stiff;
346  FloatArray alpha(6);
347  alpha.zero();
348  if ( nDimensions == 2 ) {
349  stiff = {
350  Diso->at(1, 1), Diso->at(1, 2), 0, 0, 0, Diso->at(1, 6), \
351  Diso->at(2, 2), 0, 0, 0, Diso->at(2, 6), 33, 0, 0, 0, 44, \
352  0, 0, 55, 0, Diso->at(6, 6)
353  };
354  } else if ( nDimensions == 3 ) {
355  stiff = {
356  Diso->at(1, 1), Diso->at(1, 2), Diso->at(1, 3), Diso->at(1, 4), Diso->at(1, 5), Diso->at(1, 6), \
357  Diso->at(2, 2), Diso->at(2, 3), Diso->at(2, 4), Diso->at(2, 5), Diso->at(2, 6), \
358  Diso->at(3, 3), Diso->at(3, 4), Diso->at(3, 5), Diso->at(3, 6), \
359  Diso->at(4, 4), Diso->at(4, 5), Diso->at(4, 6), \
360  Diso->at(5, 5), Diso->at(5, 6), Diso->at(6, 6)
361  };
362  } else {
363  OOFEM_ERROR("Invalid number of dimensions. Only 2d and 3d domains are supported in QC simulation. \n");
364  }
367  irMat.setField(0.0, _IFT_Material_density);
368  } else {
369  OOFEM_ERROR("Invalid homMtrxType");
370  }
371  mat->initializeFrom(& irMat);
372  d->setMaterial(nmat, mat);
373 
374 
375  // add material to interpolation elements
376  int ninterpelem = interpolationElementNumbers.giveSize();
377  int elNum;
378  int matNum = nmat;
379  // over interpolation elements
380  for ( int i = 1; i <= ninterpelem; i++ ) {
381  elNum = interpolationElementNumbers.at(i);
382  d->giveElement(elNum)->setMaterial(matNum);
383  }
384 }
385 
386 
387 void
389 // Approach 3 --- Local homogenization (for each element individually)
390 {
392 
393  // decision which links will solved explicitly/homogenized
395  elemList.zero();
397  nodeList.zero();
398 
399  bool stiffnesAssigned = false;
400  FloatMatrix stfTensorOf1Link(9, 9);
401  double s0Of1Link;
402 
403  int noIntEl = interpolationElementNumbers.giveSize();
404  std :: vector< FloatMatrix * >individualStiffnessMatrices;
405  std :: vector< FloatMatrix * >individualStiffnessTensors;
406  individualStiffnessMatrices.resize(noIntEl);
407  individualStiffnessTensors.resize(noIntEl);
408  for ( int i = 0; i <= noIntEl - 1; i++ ) {
409  individualStiffnessTensors [ i ] = new FloatMatrix(9, 9);
410  individualStiffnessTensors [ i ]->zero();
411  individualStiffnessMatrices [ i ] = new FloatMatrix(6, 6);
412  individualStiffnessMatrices [ i ]->zero();
413  }
414 
415  FloatArray individualS0;
416  individualS0.resize(noIntEl);
417  individualS0.zero();
418 
419 
420  // ovel all elements // TO DO improve by using list of links only
421  for ( int i = 1; i <= d->giveNumberOfElements(); i++ ) {
422  if ( d->giveElement(i)->giveNumberOfNodes() > 2 ) {
423  elemList.at(i) = 1; // all interpolation elements will be active
424  continue; // skip interpolation elements
425  }
426 
427  Element *e = d->giveElement(i);
428  int n1 = e->giveNode(1)->giveQcNodeType();
429  int n2 = e->giveNode(2)->giveQcNodeType();
430  // repnode--repnode - save link to elemList
431  if ( ( n1 == 1 ) && ( n2 == 1 ) ) {
432  elemList.at(i) = 1;
433  nodeList.at( e->giveNode(1)->giveNumber() ) = 1;
434  nodeList.at( e->giveNode(2)->giveNumber() ) = 1;
435 
436  // hanging--hanging
437  } else if ( ( n1 == 2 ) && ( n2 == 2 ) ) {
438  // one or both ends ale not located outside interpolation elements
439  // link will be solved explicitly
440  qcNode *qn1 = static_cast< qcNode * >( e->giveNode(1) );
441  qcNode *qn2 = static_cast< qcNode * >( e->giveNode(2) );
442  int mn = e->giveNode(1)->giveNumber();
443  int hn = e->giveNode(2)->giveNumber();
444  if ( qn1->giveMasterElementNumber() < 0 || qn2->giveMasterElementNumber() < 0 ) {
445  elemList.at(i) = 1;
446  nodeList.at(mn) = 1;
447  nodeList.at(hn) = 1;
448  }
449  // both end are inside interp. elements
450  // link will be homogenized
451  else {
452  stiffnesAssigned = false;
453  stiffnesAssigned = stiffnessAssignment(individualStiffnessTensors, individualS0, d, e, qn1, qn2);
454  // stiffnesAssigned = stiffnessAssignment( individualStiffnessTensors, individualS0 );
455  if ( !stiffnesAssigned ) {
456  // stiffness of this link was no successfully asigned to interpolation elements - link will be solved explicitly
457  elemList.at(i) = 1;
458  nodeList.at(mn) = 1;
459  nodeList.at(hn) = 1;
460  }
461  }
462 
463  // repnode--hanging
464  } else if ( ( n1 == 1 ) && ( n2 == 2 ) ) {
465  int mn = e->giveNode(1)->giveNumber();
466  int hn = e->giveNode(2)->giveNumber();
467  qcNode *qn = static_cast< qcNode * >( e->giveNode(2) );
468  int masterElem = qn->giveMasterElementNumber();
469  IntArray nodesOfMasterElem = d->giveElement(masterElem)->giveDofManArray();
470  // is masterNode one of nodes of master element?
471  if ( !nodesOfMasterElem.contains(mn) ) { // master node is not vertex of hn master elements
472  elemList.at(i) = 1;
473  nodeList.at(mn) = 1;
474  nodeList.at(hn) = 1;
475  } else {
476  // stiffness of link is assigned to masterElement
477  computeStiffnessTensorOf1Link(stfTensorOf1Link, s0Of1Link, e, d);
478  individualStiffnessTensors [ interpolationElementIndices.at(masterElem) - 1 ]->add(stfTensorOf1Link);
479  individualS0 [ interpolationElementIndices.at(masterElem) - 1 ] += s0Of1Link;
480  }
481 
482  // hanging--repnode
483  } else if ( ( n1 == 2 ) && ( n2 == 1 ) ) {
484  int mn = e->giveNode(2)->giveNumber();
485  int hn = e->giveNode(1)->giveNumber();
486  qcNode *n = static_cast< qcNode * >( e->giveNode(1) );
487  int masterElem = n->giveMasterElementNumber();
488  IntArray nodesOfMasterElem = d->giveElement(masterElem)->giveDofManArray();
489  // is masterNode one of nodes of master element?
490  if ( !nodesOfMasterElem.contains(mn) ) { // master node is not vertex of hn master elements
491  elemList.at(i) = 1;
492  nodeList.at(hn) = 1;
493  nodeList.at(mn) = 1;
494  } else {
495  // stiffness of link is assigned to masterElement
496  computeStiffnessTensorOf1Link(stfTensorOf1Link, s0Of1Link, e, d);
497  individualStiffnessTensors [ interpolationElementIndices.at(masterElem) - 1 ]->add(stfTensorOf1Link);
498  individualS0 [ interpolationElementIndices.at(masterElem) - 1 ] += s0Of1Link;
499  }
500 
501 
502  // no QcNode
503  } else {
504  OOFEM_WARNING("Element %d with non-qcNode can not be homogenized", i);
505  elemList.at(i) = 1;
506  nodeList.at( e->giveNode(1)->giveNumber() ) = 1;
507  nodeList.at( e->giveNode(2)->giveNumber() ) = 1;
508  }
509  } // ovel all links
510 
511 
512  // convert stiffness tensors to stiff matrices and divide by volume of element
513  double volumeofEl, th;
514  for ( int i = 1; i <= noIntEl; i++ ) {
515  transformStiffnessTensorToMatrix(individualStiffnessMatrices [ i - 1 ], individualStiffnessTensors [ i - 1 ]);
516  volumeofEl = d->giveElement( interpolationElementNumbers.at(i) )->computeVolumeAreaOrLength();
517  if ( d->giveNumberOfSpatialDimensions() == 2 ) {
518  th = d->giveElement( interpolationElementNumbers.at(i) )->giveCrossSection()->give(CS_Thickness, NULL);
519  volumeofEl = volumeofEl / th;
520  }
521  individualStiffnessMatrices [ i - 1 ]->times(1 / volumeofEl);
522  }
523  // ...
524 
525 
526 
527  // deactivate homogenized nodes
528  for ( int i = 1; i <= d->giveNumberOfDofManagers(); i++ ) {
529  if ( !nodeList(i) ) {
531  }
532  }
533 
534  // deactivate homogenized elements
535  for ( int i = 1; i <= d->giveNumberOfElements(); i++ ) {
536  if ( !elemList(i) ) {
538  }
539  }
540 
541  // set up lists of active nodes and elements
542  QClinearStatic *em = dynamic_cast< QClinearStatic * >( d->giveEngngModel() );
543  if ( em ) {
546  } else {
547  OOFEM_ERROR("Quasicontinuum can be applied only in QClinearStatic engngm");
548  }
549 
550 
551 
552 
553  // new materials with homogenized stiffness
554  int originalNumberOfMaterials = d->giveNumberOfMaterialModels();
555  d->resizeMaterials(originalNumberOfMaterials + noIntEl);
556  int nmat = originalNumberOfMaterials;
557  Material *mat;
558  DynamicInputRecord irMat;
559 
560  // isotropic
561  if ( homMtrxType == 1 ) {
562  double homogenizedE, homogenizedNu;
563  for ( int i = 1; i <= noIntEl; i++ ) {
564  nmat++;
565  homogenizationOfStiffMatrix(homogenizedE, homogenizedNu, individualStiffnessMatrices [ i - 1 ]);
566  mat = classFactory.createMaterial("IsoLE", nmat, d);
567  irMat.setField(homogenizedE, _IFT_IsotropicLinearElasticMaterial_e);
568  irMat.setField(homogenizedNu, _IFT_IsotropicLinearElasticMaterial_n);
570  irMat.setField(0.0, _IFT_Material_density);
571 
572  mat->initializeFrom(& irMat);
573  d->setMaterial(nmat, mat);
574  }
575 
576  // anisotropic
577  } else if ( homMtrxType == 2 ) {
578  FloatMatrix Da;
579  for ( int i = 1; i <= noIntEl; i++ ) {
580  Da = * individualStiffnessMatrices [ i - 1 ];
581  nmat++;
582  mat = classFactory.createMaterial("AnisoLE", nmat, d);
583  FloatArray stiff;
584  FloatArray alpha(6);
585  alpha.zero();
586  if ( nDimensions == 2 ) {
587  stiff = {
588  Da.at(1, 1), Da.at(1, 2), 0, 0, 0, Da.at(1, 6), \
589  Da.at(2, 2), 0, 0, 0, Da.at(2, 6), 33, 0, 0, 0, 44, \
590  0, 0, 55, 0, Da.at(6, 6)
591  };
592  } else if ( nDimensions == 3 ) {
593  stiff = {
594  Da.at(1, 1), Da.at(1, 2), Da.at(1, 3), Da.at(1, 4), Da.at(1, 5), Da.at(1, 6), \
595  Da.at(2, 2), Da.at(2, 3), Da.at(2, 4), Da.at(2, 5), Da.at(2, 6), \
596  Da.at(3, 3), Da.at(3, 4), Da.at(3, 5), Da.at(3, 6), \
597  Da.at(4, 4), Da.at(4, 5), Da.at(4, 6), \
598  Da.at(5, 5), Da.at(5, 6), Da.at(6, 6)
599  };
600  } else {
601  OOFEM_ERROR("Invalid number of dimensions. Only 2d and 3d domains are supported in QC simulation. \n");
602  }
605  irMat.setField(0.0, _IFT_Material_density);
606 
607  mat->initializeFrom(& irMat);
608  d->setMaterial(nmat, mat);
609  }
610  } else {
611  OOFEM_ERROR("Invalid homMtrxType");
612  }
613 
614 
615  // add material to interpolation elements
616  int elNum;
617  int matNum = originalNumberOfMaterials;
618  // over interpolation elements
619  for ( int i = 1; i <= noIntEl; i++ ) {
620  matNum++;
621  elNum = interpolationElementNumbers.at(i);
622  d->giveElement(elNum)->setMaterial(matNum);
623  }
624 
625 
626 
627  // delete
628  for ( unsigned int i = 0; i <= individualStiffnessTensors.size() - 1; i++ ) {
629  delete individualStiffnessTensors [ i ];
630  delete individualStiffnessMatrices [ i ];
631  }
632 
633  for ( unsigned int i = 0; i <= connectivityTable.size() - 1; i++ ) {
634  delete connectivityTable [ i ];
635  }
636 }
637 
638 
639 
640 void
641 Quasicontinuum :: homogenizationOfStiffMatrix(double &homogenizedE, double &homogenizedNu, FloatMatrix *Diso)
642 // identification of E and Nu from isotropic stiffness matrix
643 {
644  // for 2d (planme stress only)
645  if ( this->nDimensions == 2 ) {
646  double a = Diso->at(1, 1);
647  double b = Diso->at(1, 2);
648  double d = Diso->at(2, 2);
649  double f = Diso->at(6, 6); // xy stiff is in last column
650  // Norma A - using eigen vectors - not implemented in 3D
651  if ( 33 * a + 2 * b + 33 * d + 4 * f == 0 ) { // zero matrix = empty interpolation element -> zero stiffness
652  homogenizedE = 1.0e-20;
653  homogenizedNu = 0;
654  } else {
655  homogenizedE = 4 * ( a + 2 * b + d ) * ( 4 * a - 8 * b + 4 * d + f ) / ( 33 * a + 2 * b + 33 * d + 4 * f );
656  homogenizedNu = ( a + 66 * b + d - 4 * f ) / ( 33 * a + 2 * b + 33 * d + 4 * f );
657  }
658  // Norma B - sum of squares of differences (with weight functions)
659  /*
660  * if (3*a+12*b+3*d==0) {// zero matrix = empty interpolation element -> zero stiffness
661  * homogenizedE = 1.0e-20;
662  * homogenizedNu = 0;
663  * } else {
664  * homogenizedE = (a+2*b+d)*(a+14*b+d)/2/(3*a+12*b+3*d);
665  * homogenizedNu = 2*(a-b+d)/(3*a+12*b+3*d);
666  * }
667  */
668 
669  //for 3d
670  } else if ( this->nDimensions == 3 ) {
671  // Norma A - using eigen vectors
672  //- not implemented in 3D
673 
674  // Norma B - sum of squares of differences (with weight functions)
675  double a = Diso->at(1, 1) + Diso->at(2, 2) + Diso->at(3, 3);
676  double b = Diso->at(4, 4) + Diso->at(5, 5) + Diso->at(6, 6);
677 
678  if ( 5 * a + 13 * b == 0 ) { // zero matrix = empty interpolation element -> zero stiffness
679  homogenizedE = 1.0e-20;
680  homogenizedNu = 0;
681  } else {
682  homogenizedE = ( a + 2 * b ) * ( a + 11 * b ) / ( 15 * a + 39 * b );
683  homogenizedNu = ( 2 * a + b ) / ( 5 * a + 13 * b );
684  }
685  } else {
686  OOFEM_ERROR("Invalid number of dimensions. Only 2d and 3d domains are supported in QC simulation. \n");
687  }
688 }
689 
690 
691 void
693 {
694  // stiffness tensor 3^4 of 1 link represented by 9x9 mtrx
695  D1.resize(9, 9);
696  D1.zero();
697 
698  S0 = 0.;
699 
700  double Etruss = 0.;
701  double Ry = 0.;
703  if ( !qcmei ) {
704  OOFEM_ERROR("Material doesn't implement the required QC interface!");
705  }
706  Etruss = qcmei->giveQcElasticParamneter();
707  Ry = qcmei->giveQcPlasticParamneter();
708 
709  double area = 0.;
710  SimpleCrossSection *cs = dynamic_cast< SimpleCrossSection * >( e->giveCrossSection() );
711  if ( cs ) {
712  area = cs->give(CS_Area, NULL);
713  } else {
714  OOFEM_ERROR( "Invalid CrossSection of link %d. simpleCS is only supported CS of links in QC simulation. \n", e->giveGlobalNumber() );
715  }
716 
717  double x1 = e->giveDofManager(1)->giveCoordinates()->at(1);
718  double y1 = e->giveDofManager(1)->giveCoordinates()->at(2);
719  double z1 = e->giveDofManager(1)->giveCoordinates()->at(3);
720  double x2 = e->giveDofManager(2)->giveCoordinates()->at(1);
721  double y2 = e->giveDofManager(2)->giveCoordinates()->at(2);
722  double z2 = e->giveDofManager(2)->giveCoordinates()->at(3);
723  double TrussLength = sqrt( ( x2 - x1 ) * ( x2 - x1 ) + ( y2 - y1 ) * ( y2 - y1 ) + ( z2 - z1 ) * ( z2 - z1 ) );
724  double n [ 3 ];
725  n [ 0 ] = ( x2 - x1 ) / TrussLength;
726  n [ 1 ] = ( y2 - y1 ) / TrussLength;
727  n [ 2 ] = ( z2 - z1 ) / TrussLength;
728  // create stiffness tensor
729  int i, j, k, l;
730  for ( i = 1; i <= 3; i++ ) {
731  for ( j = 1; j <= 3; j++ ) {
732  for ( k = 1; k <= 3; k++ ) {
733  for ( l = 1; l <= 3; l++ ) {
734  D1.at(3 * ( i - 1 ) + k, 3 * ( j - 1 ) + l) = TrussLength * Etruss * n [ i - 1 ] * n [ j - 1 ] * n [ k - 1 ] * n [ l - 1 ];
735  }
736  }
737  }
738  }
739 
740  // homogenized yeld limit
741  if ( !std :: isinf(Ry) ) {
742  S0 += TrussLength * area * Ry / 3;
743  }
744 }
745 
746 
747 
748 void
749 Quasicontinuum :: createGlobalStiffnesMatrix(FloatMatrix *D, double &S0, Domain *d, int homMtrxType, double volumeOfInterpolationMesh)
750 //
751 {
752  int i, j, k, l, t;
753  int noe = d->giveNumberOfElements();
754 
755 
756  S0 = 0;
757  //double StiffessTensor[81]; // stiffness tensor 3^4 represented by 9x9 mtrx
758  FloatMatrix StiffTens;
759  FloatMatrix *StiffnessTensor = & StiffTens;
760  StiffnessTensor->resize(9, 9);
761  StiffnessTensor->zero();
762 
763  //double D1[81]; // stiffness tensor 3^4 of 1 link represented by 9x9 mtrx
764  FloatMatrix D1(9, 9);
765  D1.zero();
766 
767  double TrussLength, x1, x2, y1, y2, z1, z2;
768  double n [ 3 ];
769  // over all elements
770  for ( t = 1; t <= noe; t++ ) {
771  // ??km?? TO DO: use list of interpolation elements instead of skipping...
772  if ( d->giveElement(t)->giveNumberOfDofManagers() > 2 ) {
773  continue; // skip interpolation elements
774  }
775 
776  // TO DO: function "computeStiffnessMatrixOf1Link" can be used here
777  double Etruss = 0.;
778  double Ry = 0.;
780  if ( !qcmei ) {
781  OOFEM_ERROR("Material doesn't implement the required QC interface!");
782  }
783  Etruss = qcmei->giveQcElasticParamneter();
784  Ry = qcmei->giveQcPlasticParamneter();
785 
786  double area = 0.;
787  SimpleCrossSection *cs = dynamic_cast< SimpleCrossSection * >( d->giveElement(t)->giveCrossSection() );
788  if ( cs ) {
789  area = cs->give(CS_Area, NULL);
790  } else {
791  OOFEM_ERROR( "Invalid CrossSection of link %d. simpleCS is only supported CS of links in QC simulation. \n", d->giveElement(t)->giveGlobalNumber() );
792  }
793 
794  x1 = d->giveElement(t)->giveDofManager(1)->giveCoordinates()->at(1);
795  y1 = d->giveElement(t)->giveDofManager(1)->giveCoordinates()->at(2);
796  z1 = d->giveElement(t)->giveDofManager(1)->giveCoordinates()->at(3);
797  x2 = d->giveElement(t)->giveDofManager(2)->giveCoordinates()->at(1);
798  y2 = d->giveElement(t)->giveDofManager(2)->giveCoordinates()->at(2);
799  z2 = d->giveElement(t)->giveDofManager(2)->giveCoordinates()->at(3);
800  TrussLength = sqrt( ( x2 - x1 ) * ( x2 - x1 ) + ( y2 - y1 ) * ( y2 - y1 ) + ( z2 - z1 ) * ( z2 - z1 ) );
801  n [ 0 ] = ( x2 - x1 ) / TrussLength;
802  n [ 1 ] = ( y2 - y1 ) / TrussLength;
803  n [ 2 ] = ( z2 - z1 ) / TrussLength;
804  // stiffness tensor
805  for ( i = 1; i <= 3; i++ ) {
806  for ( j = 1; j <= 3; j++ ) {
807  for ( k = 1; k <= 3; k++ ) {
808  for ( l = 1; l <= 3; l++ ) {
809  D1.at(3 * ( i - 1 ) + k, 3 * ( j - 1 ) + l) = TrussLength * Etruss * n [ i - 1 ] * n [ j - 1 ] * n [ k - 1 ] * n [ l - 1 ];
810  }
811  }
812  }
813  }
814 
815  StiffnessTensor->add(D1); // add stfMtrX of 1 link to global stfMtrx
816 
817  if ( !std :: isinf(Ry) ) {
818  S0 += TrussLength * area * Ry / 3;
819  }
820  } // for t=1:note
821 
822  // divide by volume of element
823  StiffnessTensor->times(1 / volumeOfInterpolationMesh);
824  S0 = S0 / volumeOfInterpolationMesh;
825 
826  transformStiffnessTensorToMatrix(D, StiffnessTensor);
827 }
828 
829 
830 bool
831 Quasicontinuum :: stiffnessAssignment(std :: vector< FloatMatrix * > &individualStiffnessTensors, FloatArray &individualS0, Domain *d, Element *e, qcNode *qn1, qcNode *qn2)
832 {
833  register int i;
834  FloatMatrix stfTensorOf1Link(9, 9), contribution(9, 9);
835  double s0Of1Link;
836  int nel, indx;
837 
838 
839  // numbers of elements containing ending nodes
840  int end1 = qn1->giveMasterElementNumber();
841  int end2 = qn2->giveMasterElementNumber();
842 
843  // Coords and length
844  FloatArray *coords1 = qn1->giveCoordinates();
845  FloatArray *coords2 = qn2->giveCoordinates();
846  //coordsDif->subtract(*coords1);
847  FloatArray coordsDif(2);
848  coordsDif.beDifferenceOf(* coords1, * coords2);
849  double totalLength = coordsDif.computeNorm();
850 
851  /*
852  * if ( (end1==0) && (end2==0) ) { // both ends are located outside
853  * // this case is solved explicitly before StiffnessAssignment
854  * } elseif ( (end1==0) || (end2==0) ) { // both ends are located outside
855  * // this case is solved explicitly before StiffnessAssignment
856  * }
857  */
858 
859  if ( end1 == end2 ) { // both ends are in the same el.
860  // stiffness of link is assigned to this el.
861  computeStiffnessTensorOf1Link(stfTensorOf1Link, s0Of1Link, e, d);
862  nel = end1; // number of interp element
863  indx = interpolationElementIndices.at(nel); // number in list of interp elem
864  individualStiffnessTensors [ indx - 1 ]->add(stfTensorOf1Link);
865  individualS0 [ indx - 1 ] += s0Of1Link;
866 
867  return true;
868  } else { // ends are located in different elements -> all intersected el needs to be found
869  IntArray intersected; // to store numbers of intersected elem
870  FloatArray lengths; // to store lengths ofintersections
871  intersected.clear();
872  lengths.clear();
873 
874  computeIntersectionsOfLinkWithInterpElements(intersected, lengths, d, e, qn1, qn2);
875 
876  // check: is there any intersectted elements
877  int noIntersected = lengths.giveSize();
878  if ( noIntersected == 0 ) {
879  return false;
880  }
881 
882  // check: is stiff. assigned to all parts of the link?
883  // TO DO: the same check is done in "computeIntersectionsOfLinkWith...".
884  // Use only return value here
885  lengths.times(1 / totalLength); // it is better to use relative lengths
886  double sumLength = 0.0;
887  for ( int i = 1; i <= noIntersected; i++ ) {
888  sumLength += lengths [ i - 1 ];
889  }
890  if ( sumLength > 1.0001 ) {
891  return false;
892  } else if ( sumLength < 0.9999 ) {
893  return false;
894  }
895 
896  // if check OK ( sumLength == 1 )
897  computeStiffnessTensorOf1Link(stfTensorOf1Link, s0Of1Link, e, d);
898 
899  for ( i = 1; i <= noIntersected; i++ ) { // add stiff to all intersected elem
900  nel = intersected [ i - 1 ]; // number of element
901  indx = interpolationElementIndices.at(nel); // number in list of interp elem
902  contribution = stfTensorOf1Link;
903  contribution.times( lengths.at(i) );
904  individualStiffnessTensors [ indx - 1 ]->add(contribution);
905  individualS0 [ indx - 1 ] += lengths.at(i) * s0Of1Link;
906  }
907  return true;
908  } // ends are located in different elements
909 
910  return false;
911 }
912 
913 
914 void
916 {
917  int dim = d->giveNumberOfSpatialDimensions();
918  if ( dim == 2 ) {
919  computeIntersectionsOfLinkWith2DTringleElements(intersected, lengths, d, e, qn1, qn2);
920  } else if ( dim == 3 ) {
921  computeIntersectionsOfLinkWith3DTetrahedraElements(intersected, lengths, d, e, qn1, qn2);
922  } else {
923  OOFEM_ERROR("unsupported number of dimensions");
924  }
925 }
926 
927 bool
929 {
930  intersected.clear();
931  lengths.clear();
932 
933  // coordinates of ending nodes
934  FloatArray *X1 = qn1->giveCoordinates();
935  FloatArray *X2 = qn2->giveCoordinates();
936 
937  double TotalLength = sqrt( ( X2->at(1) - X1->at(1) ) * ( X2->at(1) - X1->at(1) ) + ( X2->at(2) - X1->at(2) ) * ( X2->at(2) - X1->at(2) ) );
938 
939  // numbers of elements containing ending nodes
940  int end1 = qn1->giveMasterElementNumber();
941  int end2 = qn2->giveMasterElementNumber();
942 
943  if ( end2 < end1 ) { // swap?
944  X2 = qn1->giveCoordinates();
945  X1 = qn2->giveCoordinates();
946  end2 = qn1->giveMasterElementNumber();
947  end1 = qn2->giveMasterElementNumber();
948  }
949 
950  int numberOfIntersected = 0; // number of intersected elements
951 
952  FloatArray *A, *B, *C, vector(2);
953  FloatArray intersectCoordsX, intersectCoordsY;
954  double iLen, sumLength = 0.0;
955  register int nn, ii, i;
956 
957  // ftiffness assignment for starting element (end1)
958  int iel = end1;
959  FloatArray iLens;
960 
961  A = d->giveElement(iel)->giveNode(1)->giveCoordinates();
962  B = d->giveElement(iel)->giveNode(2)->giveCoordinates();
963  C = d->giveElement(iel)->giveNode(3)->giveCoordinates();
964  numberOfIntersected = intersectionTestSegmentTriangle2D(intersectCoordsX, intersectCoordsY, A, B, C, X1, X2);
965 
966  switch ( numberOfIntersected ) {
967  case 0: // no intersection with this element
968  lengths.push_back(0);
969  intersected.followedBy(iel); // cislo elementu
970  break;
971  case 1: // link starts in this element
972  vector.at(1) = intersectCoordsX.at(1) - X1->at(1);
973  vector.at(2) = intersectCoordsY.at(1) - X1->at(2);
974  iLen = vector.computeNorm();
975  if ( iLen / TotalLength >= 1.e-6 ) { // intersection is not negligible
976  lengths.push_back(iLen);
977  intersected.followedBy(iel);
978  } else { // intersection is negligible
979  lengths.push_back(0);
980  intersected.followedBy(iel); // cislo elementu
981  }
982  break;
983  default: // 2 or more intersections
984 
985  //lengths of all possible cobmination of intersection points
986  iLens.resize(numberOfIntersected);
987  for ( int nn = 1; nn <= numberOfIntersected; nn++ ) {
988  vector.at(1) = intersectCoordsX.at(nn) - X1->at(1);
989  vector.at(2) = intersectCoordsY.at(nn) - X1->at(2);
990  iLens.at(nn) = vector.computeNorm();
991  }
992  iLen = iLens.at( iLens.giveIndexMaxElem() ); // real length is maximal one
993  if ( iLen / TotalLength >= 1.e-6 ) { // intersection is not negligible
994  lengths.push_back(iLen);
995  intersected.followedBy(iel);
996  } else { // intersection is negligible
997  lengths.push_back(0);
998  intersected.followedBy(iel); // cislo elementu
999  }
1000  break;
1001  }
1002 
1003 
1004  // ftiffness assignment for the rest of the link
1005  bool breakFor = false;
1006  IntArray neighboursList, alreadyTestedElemList;
1007  alreadyTestedElemList.clear();
1008 
1010 
1011  //bool secondEndReached = false;
1012  int nextElement = end1;
1013  int nel, index, m;
1014  register int k;
1015 
1016  for ( ii = 1; ii <= nome; ii++ ) { // searching for nex intersected element
1017  nel = nextElement; // previous element number
1018  nextElement = 0; // next element number (to be found)
1019 
1020  // list of connected elements
1021  neighboursList = * connectivityTable [ interpolationElementIndices.at(nel) - 1 ];
1022 
1023  // delete already intersected elements
1024  for ( k = 1; k <= intersected.giveSize(); k++ ) {
1025  index = neighboursList.findFirstIndexOf( intersected.at(k) );
1026  if ( index != 0 ) {
1027  neighboursList.erase(index);
1028  }
1029  }
1030 
1031  // delete already tested elements
1032  for ( k = 1; k <= alreadyTestedElemList.giveSize(); k++ ) {
1033  index = neighboursList.findFirstIndexOf( alreadyTestedElemList.at(k) );
1034  if ( index != 0 ) {
1035  neighboursList.erase(index);
1036  }
1037  }
1038 
1039 
1040  if ( neighboursList.giveSize() == 0 ) { // no elements to continue
1041  //secondEndReached = false;
1042  break; // ends searching for nex intersected element
1043  }
1044 
1045  for ( k = 1; k <= neighboursList.giveSize(); k++ ) { // testing of neighbouring elements
1046  iel = neighboursList.at(k);
1047 
1048  A = d->giveElement(iel)->giveNode(1)->giveCoordinates();
1049  B = d->giveElement(iel)->giveNode(2)->giveCoordinates();
1050  C = d->giveElement(iel)->giveNode(3)->giveCoordinates();
1051  numberOfIntersected = intersectionTestSegmentTriangle2D(intersectCoordsX, intersectCoordsY, A, B, C, X1, X2);
1052 
1053  breakFor = false; // go to next element?
1054  switch ( numberOfIntersected ) {
1055  case 0: // no intersection with (iel) element -> go to next element
1056  alreadyTestedElemList.followedBy(iel);
1057  continue;
1058  case 1: // link ends here or there in only fake intersection
1059  if ( iel == end2 ) { // end of the link in in this element
1060  vector.at(1) = intersectCoordsX.at(1) - X2->at(1);
1061  vector.at(2) = intersectCoordsY.at(1) - X2->at(2);
1062  iLen = vector.computeNorm();
1063  if ( iLen / TotalLength >= 1.e-6 ) { // intersection is not negligible
1064  lengths.push_back(iLen);
1065  intersected.followedBy(iel);
1066  // no break here, some elements still needs to be tested
1067  } else { // intersection is negligible
1068  lengths.push_back(0);
1069  intersected.followedBy(iel); // cislo elementu
1070  }
1071  } else { // fake intersection (only touch with edge)
1072  lengths.push_back(0);
1073  intersected.followedBy(iel);
1074  }
1075  break;
1076  default: // 2 or more intersections
1077  if ( iel == end2 ) { // end of the link is inthis element
1078  //lengths of all possible cobmination of intersection points
1079  iLens.resize(numberOfIntersected);
1080  for ( int nn = 1; nn <= numberOfIntersected; nn++ ) {
1081  vector.at(1) = intersectCoordsX.at(nn) - X2->at(1);
1082  vector.at(2) = intersectCoordsY.at(nn) - X2->at(2);
1083  iLens.at(nn) = vector.computeNorm();
1084  }
1085  iLen = iLens.at( iLens.giveIndexMaxElem() ); // real length is maximal
1086  if ( iLen / TotalLength >= 1.e-6 ) { // intersection is not negligible
1087  lengths.push_back(iLen);
1088  intersected.followedBy(iel);
1089  nextElement = iel;
1090  breakFor = true;
1091  break; // go to next element
1092  } else { // intersection is negligible
1093  lengths.push_back(0);
1094  intersected.followedBy(iel);
1095  }
1096  } else { // end of the link is NOT inthis elem
1097  iLens.resize(3); // combinations: 12 13 23
1098  iLens.zero();
1099  m = 0;
1100  for ( nn = 1; nn < numberOfIntersected; nn++ ) {
1101  for ( ii = nn + 1; ii <= numberOfIntersected; ii++ ) {
1102  m += 1;
1103  vector.at(1) = intersectCoordsX.at(nn) - intersectCoordsX.at(ii);
1104  vector.at(2) = intersectCoordsY.at(nn) - intersectCoordsY.at(ii);
1105  iLens.at(m) = vector.computeNorm();
1106  }
1107  }
1108  iLen = iLens.at( iLens.giveIndexMaxElem() ); // real length is maximum
1109  if ( iLen / TotalLength >= 1.e-6 ) { // intersection is not negligible
1110  lengths.push_back(iLen);
1111  intersected.followedBy(iel);
1112  nextElement = iel;
1113  breakFor = true;
1114  break; // go to next element
1115  } else { // intersection is negligible
1116  lengths.push_back(0);
1117  intersected.followedBy(iel);
1118  }
1119  } // end of the link is inthis element
1120  } // switch
1121 
1122  if ( breakFor ) {
1123  break;
1124  } // go to next interpolation element
1125  } // over neigbouring elements
1126 
1127 
1128  // check: all length of link is assembles
1129  sumLength = 0.0;
1130  for ( i = 1; i <= lengths.giveSize(); i++ ) {
1131  sumLength += lengths.at(i);
1132  }
1133  if ( ( fabs(sumLength / TotalLength - 1) ) <= ( 0.0001 ) ) {
1134  return true;
1135  }
1136 
1137  // check: there is no neighnouring element to continue
1138  if ( nextElement == 0 ) {
1139  return false;
1140  }
1141  } // searching for nex intersected element
1142 
1143  return false;
1144 }
1145 
1146 
1147 
1148 
1149 
1150 bool
1152 {
1153  intersected.clear();
1154  lengths.clear();
1155 
1156  // coordinates of ending nodes
1157  FloatArray *X1 = qn1->giveCoordinates();
1158  FloatArray *X2 = qn2->giveCoordinates();
1159 
1160  double TotalLength = sqrt( ( X2->at(1) - X1->at(1) ) * ( X2->at(1) - X1->at(1) ) + ( X2->at(2) - X1->at(2) ) * ( X2->at(2) - X1->at(2) ) + ( X2->at(3) - X1->at(3) ) * ( X2->at(3) - X1->at(3) ) );
1161 
1162  // numbers of elements containing ending nodes
1163  int end1 = qn1->giveMasterElementNumber();
1164  int end2 = qn2->giveMasterElementNumber();
1165 
1166  if ( end2 < end1 ) { // swap?
1167  X2 = qn1->giveCoordinates();
1168  X1 = qn2->giveCoordinates();
1169  end2 = qn1->giveMasterElementNumber();
1170  end1 = qn2->giveMasterElementNumber();
1171  }
1172 
1173  int numberOfIntersected = 0; // number of intersected elements
1174 
1175  FloatArray *A, *B, *C, *D, vector(3);
1176  FloatArray intersectCoordsX, intersectCoordsY, intersectCoordsZ;
1177  double iLen, sumLength = 0.0;
1178  register int nn, ii, i;
1179 
1180 
1181 
1182  // ftiffness assignment for starting element (end1)
1183 
1184  int iel = end1;
1185  FloatArray iLens;
1186 
1187  A = d->giveElement(iel)->giveNode(1)->giveCoordinates();
1188  B = d->giveElement(iel)->giveNode(2)->giveCoordinates();
1189  C = d->giveElement(iel)->giveNode(3)->giveCoordinates();
1190  D = d->giveElement(iel)->giveNode(4)->giveCoordinates();
1191  numberOfIntersected = intersectionTestSegmentTetrahedra3D(intersectCoordsX, intersectCoordsY, intersectCoordsZ, A, B, C, D, X1, X2);
1192 
1193  switch ( numberOfIntersected ) {
1194  case 0: // no intersection with this element
1195  lengths.push_back(0);
1196  intersected.followedBy(iel); // cislo elementu
1197  break;
1198  case 1: // link starts in this element
1199  vector.at(1) = intersectCoordsX.at(1) - X1->at(1);
1200  vector.at(2) = intersectCoordsY.at(1) - X1->at(2);
1201  vector.at(3) = intersectCoordsZ.at(1) - X1->at(3);
1202  iLen = vector.computeNorm();
1203  if ( iLen / TotalLength >= 1.e-6 ) { // intersection is not negligible
1204  lengths.push_back(iLen);
1205  intersected.followedBy(iel);
1206  } else { // intersection is negligible
1207  lengths.push_back(0);
1208  intersected.followedBy(iel);
1209  }
1210  break;
1211  default: // 2 or more intersections
1212 
1213  //lengths of all possible cobmination of intersection points
1214  iLens.resize(numberOfIntersected);
1215  for ( int nn = 1; nn <= numberOfIntersected; nn++ ) {
1216  vector.at(1) = intersectCoordsX.at(nn) - X1->at(1);
1217  vector.at(2) = intersectCoordsY.at(nn) - X1->at(2);
1218  vector.at(3) = intersectCoordsZ.at(nn) - X1->at(3);
1219  iLens.at(nn) = vector.computeNorm();
1220  }
1221  iLen = iLens.at( iLens.giveIndexMaxElem() ); // real length is maximal one
1222  if ( iLen / TotalLength >= 1.e-6 ) { // intersection is not negligible
1223  lengths.push_back(iLen);
1224  intersected.followedBy(iel);
1225  } else { // intersection is negligible
1226  lengths.push_back(0);
1227  intersected.followedBy(iel);
1228  }
1229  break;
1230  }
1231 
1232 
1233 
1234  // ftiffness assignment for the rest of the link
1235 
1236  bool breakFor = false;
1237  IntArray neighboursList, alreadyTestedElemList;
1238  alreadyTestedElemList.clear();
1239 
1241 
1242  //bool secondEndReached = false;
1243  int nextElement = end1;
1244  int nel, index, m;
1245  register int k;
1246 
1247  for ( ii = 1; ii <= nome; ii++ ) { // searching for nex intersected element
1248  nel = nextElement; // previous element number
1249  nextElement = 0; // next element number (to be found)
1250 
1251  // list of connected elements
1252  neighboursList = * connectivityTable [ interpolationElementIndices.at(nel) - 1 ];
1253 
1254  // delete already intersected elements
1255  for ( k = 1; k <= intersected.giveSize(); k++ ) {
1256  index = neighboursList.findFirstIndexOf( intersected.at(k) );
1257  if ( index != 0 ) {
1258  neighboursList.erase(index);
1259  }
1260  }
1261 
1262  // delete already tested elements
1263  for ( k = 1; k <= alreadyTestedElemList.giveSize(); k++ ) {
1264  index = neighboursList.findFirstIndexOf( alreadyTestedElemList.at(k) );
1265  if ( index != 0 ) {
1266  neighboursList.erase(index);
1267  }
1268  }
1269 
1270 
1271  if ( neighboursList.giveSize() == 0 ) { // no elements to continue
1272  //secondEndReached = false;
1273  break; // ends searching for nex intersected element
1274  }
1275 
1276  for ( k = 1; k <= neighboursList.giveSize(); k++ ) { // testing of neighbouring elements
1277  iel = neighboursList.at(k);
1278 
1279  A = d->giveElement(iel)->giveNode(1)->giveCoordinates();
1280  B = d->giveElement(iel)->giveNode(2)->giveCoordinates();
1281  C = d->giveElement(iel)->giveNode(3)->giveCoordinates();
1282  D = d->giveElement(iel)->giveNode(4)->giveCoordinates();
1283  numberOfIntersected = intersectionTestSegmentTetrahedra3D(intersectCoordsX, intersectCoordsY, intersectCoordsZ, A, B, C, D, X1, X2);
1284 
1285 
1286  breakFor = false; // go to next interp elem
1287  switch ( numberOfIntersected ) {
1288  case 0: // no intersection with (iel) element -> go to next element
1289  alreadyTestedElemList.followedBy(iel);
1290  continue; // zadna tuhost se prvku neprirazuje
1291  case 1: // link ends here or there in only fake intersection
1292  if ( iel == end2 ) { // end of the link in in this element
1293  vector.at(1) = intersectCoordsX.at(1) - X2->at(1);
1294  vector.at(2) = intersectCoordsY.at(1) - X2->at(2);
1295  vector.at(3) = intersectCoordsZ.at(1) - X2->at(3);
1296  iLen = vector.computeNorm();
1297  if ( iLen / TotalLength >= 1.e-6 ) { // intersection is not negligible
1298  lengths.push_back(iLen);
1299  intersected.followedBy(iel);
1300  // no break here, some elements still needs to be tested
1301  } else { // intersection is negligible
1302  lengths.push_back(0);
1303  intersected.followedBy(iel); // cislo elementu
1304  }
1305  } else { // fake intersection (only touch with edge)
1306  lengths.push_back(0);
1307  intersected.followedBy(iel);
1308  }
1309  break;
1310  default: // 2 or more intersections
1311  if ( iel == end2 ) { // end of the link is inthis element
1312  //lengths of all possible cobmination of intersection points
1313  iLens.resize(numberOfIntersected);
1314  for ( int nn = 1; nn <= numberOfIntersected; nn++ ) {
1315  vector.at(1) = intersectCoordsX.at(nn) - X2->at(1);
1316  vector.at(2) = intersectCoordsY.at(nn) - X2->at(2);
1317  vector.at(3) = intersectCoordsZ.at(nn) - X2->at(3);
1318  iLens.at(nn) = vector.computeNorm();
1319  }
1320  iLen = iLens.at( iLens.giveIndexMaxElem() ); // real length is maximal
1321  if ( iLen / TotalLength >= 1.e-6 ) { // intersection is not negligible
1322  lengths.push_back(iLen);
1323  intersected.followedBy(iel);
1324  nextElement = iel;
1325  breakFor = true;
1326  break; // go to next element
1327  } else { // intersection is negligible
1328  lengths.push_back(0);
1329  intersected.followedBy(iel);
1330  }
1331  } else { // end of the link is NOT inthis elem
1332  iLens.resize(6); // combinations: 12 13 14 23 24 34
1333  iLens.zero();
1334  m = 0;
1335  for ( nn = 1; nn <= numberOfIntersected; nn++ ) {
1336  for ( ii = nn + 1; ii <= numberOfIntersected; ii++ ) {
1337  m += 1;
1338  vector.at(1) = intersectCoordsX.at(nn) - intersectCoordsX.at(ii);
1339  vector.at(2) = intersectCoordsY.at(nn) - intersectCoordsY.at(ii);
1340  vector.at(3) = intersectCoordsZ.at(nn) - intersectCoordsZ.at(ii);
1341  iLens.at(m) = vector.computeNorm();
1342  }
1343  }
1344  iLen = iLens.at( iLens.giveIndexMaxElem() ); // real length is maximum
1345  if ( iLen / TotalLength >= 1.e-6 ) { // intersection is not negligible
1346  lengths.push_back(iLen);
1347  intersected.followedBy(iel);
1348  nextElement = iel;
1349  breakFor = true;
1350  break; // go to next element
1351  } else { // intersection is negligible
1352  lengths.push_back(0);
1353  intersected.followedBy(iel);
1354  }
1355  } // end of the link is inthis element
1356  } // switch
1357 
1358  if ( breakFor ) {
1359  break;
1360  } // go to next interp elem
1361  } // over neigbouring elements
1362 
1363 
1364  // check: all length of link is assembles
1365  sumLength = 0.0;
1366  for ( i = 1; i <= lengths.giveSize(); i++ ) {
1367  sumLength += lengths.at(i);
1368  }
1369  if ( ( fabs(sumLength / TotalLength - 1) ) <= ( 0.0001 ) ) {
1370  return true;
1371  }
1372 
1373  // check: there is no neighnouring element to continue
1374  if ( nextElement == 0 ) {
1375  return false;
1376  }
1377  } // searching for nex intersected element
1378 
1379  return false;
1380 }
1381 
1382 void
1384 {
1385  // for each element save elements with shared node (edge, face,,,)
1386  register int i, j, k, l;
1388  int elNum1, elNum2, noVertices1, noVertices2;
1389  Element *e1, *e2;
1390  IntArray vertices1;
1391  connectivityTable.clear();
1392  connectivityTable.resize(nome);
1393 
1394 
1395 
1396  for ( j = 1; j <= nome; j++ ) { // all interpolation elements
1397  connectivityTable [ j - 1 ] = new IntArray();
1398  elNum1 = interpolationElementNumbers.at(j);
1399  e1 = d->giveElement(elNum1);
1400 
1401  // vector of vertices numbers of el j
1402  noVertices1 = e1->giveNumberOfNodes();
1403  vertices1.resize(noVertices1);
1404  for ( k = 1; k <= noVertices1; k++ ) { // fill in vector
1405  vertices1.at(k) = e1->giveNode(k)->giveNumber();
1406  }
1407  for ( i = 1; i <= nome; i++ ) {
1408  if ( i == j ) {
1409  continue;
1410  } // is not itself neighbour
1411  // vector of vertices numbers of el j
1412  elNum2 = interpolationElementNumbers.at(i);
1413  e2 = d->giveElement(elNum2);
1414  noVertices2 = e2->giveNumberOfNodes();
1415  for ( l = 1; l <= noVertices2; l++ ) { // fill in vector
1416  if ( vertices1.contains( e2->giveNode(l)->giveNumber() ) ) {
1417  connectivityTable [ j - 1 ]->followedBy(elNum2); // add number of neighbour
1418  break;
1419  }
1420  }
1421  }
1422  }
1423 }
1424 
1425 
1426 
1427 
1428 bool
1430 {
1431  // test if line (X2-X1) intersects tringle ABC (for 3D!)
1432  intersectCoords.clear();
1433 
1434  // directional vektor
1435  double L = sqrt( ( X2->at(1) - X1->at(1) ) * ( X2->at(1) - X1->at(1) ) + ( X2->at(2) - X1->at(2) ) * ( X2->at(2) - X1->at(2) ) + ( X2->at(3) - X1->at(3) ) * ( X2->at(3) - X1->at(3) ) );
1436  double Ur1 = ( X2->at(1) - X1->at(1) ) / L;
1437  double Ur2 = ( X2->at(2) - X1->at(2) ) / L;
1438  double Ur3 = ( X2->at(3) - X1->at(3) ) / L;
1439  // cross product (needed for Plucker. coords)
1440  double Vr1 = Ur2 * X1->at(3) - Ur3 *X1->at(2);
1441  double Vr2 = Ur3 * X1->at(1) - Ur1 *X1->at(3);
1442  double Vr3 = Ur1 * X1->at(2) - Ur2 *X1->at(1);
1443  // Plucker. coords of line {U,UxX1}
1444  // Pr = [Ur, Vr];
1445 
1446  // vectors of edges: (order is important!)
1447  FloatArray U1(3), U2(3), U3(3);
1448  for ( int i = 1; i <= 3; i++ ) {
1449  U1.at(i) = B->at(i) - A->at(i);
1450  U2.at(i) = C->at(i) - B->at(i);
1451  U3.at(i) = A->at(i) - C->at(i);
1452  }
1453 
1454  // cross product (needed for Plucker. coords)
1455  FloatArray V1(3), V2(3), V3(3);
1456  // cross(U1,A);
1457  V1.at(1) = U1.at(2) * A->at(3) - U1.at(3) * A->at(2);
1458  V1.at(2) = U1.at(3) * A->at(1) - U1.at(1) * A->at(3);
1459  V1.at(3) = U1.at(1) * A->at(2) - U1.at(2) * A->at(1);
1460  // cross(U2,B);
1461  V2.at(1) = U2.at(2) * B->at(3) - U2.at(3) * B->at(2);
1462  V2.at(2) = U2.at(3) * B->at(1) - U2.at(1) * B->at(3);
1463  V2.at(3) = U2.at(1) * B->at(2) - U2.at(2) * B->at(1);
1464  // cross(U3,C);
1465  V3.at(1) = U3.at(2) * C->at(3) - U3.at(3) * C->at(2);
1466  V3.at(2) = U3.at(3) * C->at(1) - U3.at(1) * C->at(3);
1467  V3.at(3) = U3.at(1) * C->at(2) - U3.at(2) * C->at(1);
1468 
1469 
1470  // permuted inner products Pr(.)P1 = Ur.V1 + U1.Vr
1471  double S1 = Ur1 * V1.at(1) + Ur2 *V1.at(2) + Ur3 *V1.at(3) + U1.at(1) * Vr1 + U1.at(2) * Vr2 + U1.at(3) * Vr3; // dot(Ur,V1) + dot(U1,Vr);
1472  double S2 = Ur1 * V2.at(1) + Ur2 *V2.at(2) + Ur3 *V2.at(3) + U2.at(1) * Vr1 + U2.at(2) * Vr2 + U2.at(3) * Vr3; // dot(Ur,V2) + dot(U2,Vr);
1473  double S3 = Ur1 * V3.at(1) + Ur2 *V3.at(2) + Ur3 *V3.at(3) + U3.at(1) * Vr1 + U3.at(2) * Vr2 + U3.at(3) * Vr3; // dot(Ur,V3) + dot(U3,Vr);
1474 
1475  // tolerance ??
1476  double TOL = 1e-11;
1477  if ( fabs(S1) <= TOL ) {
1478  S1 = 0;
1479  }
1480  if ( fabs(S2) <= TOL ) {
1481  S2 = 0;
1482  }
1483  if ( fabs(S3) <= TOL ) {
1484  S3 = 0;
1485  }
1486 
1487  // in planar test
1488  if ( ( S1 == 0 ) && ( S2 == 0 ) && ( S3 == 0 ) ) {
1489  // segment is in the same plane as triangle
1490  // in this case it means no ontersection
1491  // for tetrahedra elements there will be another 2 intersections with remaining faces
1492  intersectCoords.clear();
1493  return false; // no intersection
1494  }
1495 
1496 
1497  // same orientation -> intersection with line (not necessarily with segment)
1498  if ( ( ( S1 >= 0 ) && ( S2 >= 0 ) && ( S3 >= 0 ) ) || ( ( S1 <= 0 ) && ( S2 <= 0 ) && ( S3 <= 0 ) ) ) {
1499  // calculation of intersection coords
1500  double S = S1 + S2 + S3;
1501  // from barycentric coords
1502  double u1 = S1 / S;
1503  double u2 = S2 / S;
1504  double u3 = S3 / S;
1505  // to cartesian
1506  // Px = [x y z] = u1*C + u2*A +u3*B; // multiple by opposite vertex
1507  double x = u1 * C->at(1) + u2 *A->at(1) + u3 *B->at(1);
1508  double y = u1 * C->at(2) + u2 *A->at(2) + u3 *B->at(2);
1509  double z = u1 * C->at(3) + u2 *A->at(3) + u3 *B->at(3);
1510 
1511  // is interection on line segment
1512 
1513  // intersection parameter
1514  double t;
1515  // Pk = [x y z] = X1 + t*(X2-X1)
1516  if ( ( X2->at(1) - X1->at(1) ) != 0 ) { // x-coord will be used
1517  t = ( x - X1->at(1) ) / ( X2->at(1) - X1->at(1) ); // parameter on segment
1518  } else {
1519  if ( ( X2->at(2) - X1->at(2) ) != 0 ) { // y-coord will be used
1520  t = ( y - X1->at(2) ) / ( X2->at(2) - X1->at(2) ); // parameter on segment
1521  } else { // z-coord will be used
1522  t = ( z - X1->at(3) ) / ( X2->at(3) - X1->at(3) ); // parameter on segment
1523  }
1524  }
1525 
1526  // is interection on line segment? (0<t<1)
1527  double EPS = 1.0e-12; // numerical tolerance
1528  if ( ( 0 < t + EPS ) && ( t - EPS < 1 ) ) {
1529  intersectCoords.resize(3);
1530  intersectCoords.at(1) = x;
1531  intersectCoords.at(2) = y;
1532  intersectCoords.at(3) = z;
1533  return true;
1534  } else {
1535  intersectCoords.clear();
1536  return false; // no intersection
1537  }
1538  } else { //different orientation -> no intersection
1539  intersectCoords.clear();
1540  return false; // no intersection
1541  } // if same orientation
1542 }
1543 
1544 
1545 
1546 
1547 int
1549  // test if line (X2-X1) intersects tetrahedra ABCD (for 3D!)
1550 
1551  bool test;
1552  int numberOfIntersections = 0;
1553  intersectCoordsX.clear();
1554  intersectCoordsY.clear();
1555  intersectCoordsZ.clear();
1556  FloatArray intersectCoords;
1557 
1558  // face 1
1559  test = intersectionTestSegmentTrianglePlucker3D(intersectCoords, A, B, C, X1, X2);
1560  if ( test ) {
1561  numberOfIntersections += 1;
1562  intersectCoordsX.push_back( intersectCoords.at(1) );
1563  intersectCoordsY.push_back( intersectCoords.at(2) );
1564  intersectCoordsZ.push_back( intersectCoords.at(3) );
1565  }
1566 
1567  // face 2
1568  test = intersectionTestSegmentTrianglePlucker3D(intersectCoords, A, D, B, X1, X2);
1569  if ( test ) {
1570  numberOfIntersections += 1;
1571  intersectCoordsX.push_back( intersectCoords.at(1) );
1572  intersectCoordsY.push_back( intersectCoords.at(2) );
1573  intersectCoordsZ.push_back( intersectCoords.at(3) );
1574  }
1575 
1576  // face 3
1577  test = intersectionTestSegmentTrianglePlucker3D(intersectCoords, B, D, C, X1, X2);
1578  if ( test ) {
1579  numberOfIntersections += 1;
1580  intersectCoordsX.push_back( intersectCoords.at(1) );
1581  intersectCoordsY.push_back( intersectCoords.at(2) );
1582  intersectCoordsZ.push_back( intersectCoords.at(3) );
1583  }
1584 
1585  // face 4
1586  test = intersectionTestSegmentTrianglePlucker3D(intersectCoords, A, C, D, X1, X2);
1587  if ( test ) {
1588  numberOfIntersections += 1;
1589  intersectCoordsX.push_back( intersectCoords.at(1) );
1590  intersectCoordsY.push_back( intersectCoords.at(2) );
1591  intersectCoordsZ.push_back( intersectCoords.at(3) );
1592  }
1593 
1594  return numberOfIntersections;
1595 }
1596 
1597 
1598 
1599 int
1601  // test if line (X2-X1) intersects tringle ABC (for 2D!)
1602  bool test;
1603  int numberOfIntersections = 0;
1604  intersectCoordsX.clear();
1605  intersectCoordsY.clear();
1606  FloatArray intersectCoords;
1607 
1608  // edge 1
1609  test = intersectionTestSegmentSegment2D(intersectCoords, A, B, U1, U2);
1610  if ( test ) {
1611  numberOfIntersections += 1;
1612  intersectCoordsX.push_back( intersectCoords.at(1) );
1613  intersectCoordsY.push_back( intersectCoords.at(2) );
1614  }
1615 
1616  // edge 2
1617  test = intersectionTestSegmentSegment2D(intersectCoords, B, C, U1, U2);
1618  if ( test ) {
1619  numberOfIntersections += 1;
1620  intersectCoordsX.push_back( intersectCoords.at(1) );
1621  intersectCoordsY.push_back( intersectCoords.at(2) );
1622  }
1623 
1624  // edge 3
1625  test = intersectionTestSegmentSegment2D(intersectCoords, A, C, U1, U2);
1626  if ( test ) {
1627  numberOfIntersections += 1;
1628  intersectCoordsX.push_back( intersectCoords.at(1) );
1629  intersectCoordsY.push_back( intersectCoords.at(2) );
1630  }
1631 
1632  return numberOfIntersections;
1633 }
1634 
1635 
1636 bool
1638  //intersectionTest of 2 segments A1-A2 B1-B2 (2D only)
1639 
1640  double a1x = A1->at(1);
1641  double a1y = A1->at(2);
1642  double a2x = A2->at(1);
1643  double a2y = A2->at(2);
1644  double b1x = B1->at(1);
1645  double b1y = B1->at(2);
1646  double b2x = B2->at(1);
1647  double b2y = B2->at(2);
1648 
1649  // y=f*x+g
1650  double f1 = std :: numeric_limits< float > :: infinity();
1651  double f2 = std :: numeric_limits< float > :: infinity();
1652  if ( a1x != a2x ) {
1653  f1 = -( a2y - a1y ) / ( a1x - a2x );
1654  }
1655  if ( b1x != b2x ) {
1656  f2 = -( b2y - b1y ) / ( b1x - b2x );
1657  }
1658 
1659  double g1 = a1y - f1 * a1x;
1660  double g2 = b1y - f2 * b1x;
1661 
1662  // % parallels
1663  if ( f1 == f2 ) {
1664  intersectCoords.clear();
1665  return false;
1666  }
1667 
1668  // intersection...
1669  double x, y;
1670  if ( a1x == a2x ) {
1671  x = a1x;
1672  y = f2 * a1x + g2;
1673  } else if ( b1x == b2x ) {
1674  x = b1x;
1675  y = f1 * b1x + g1;
1676  } else {
1677  x = ( g2 - g1 ) / ( f1 - f2 );
1678  y = f1 * x + g1;
1679  }
1680 
1681  // is intersection on both segments? (it is not enough to test x-coords only)
1682  double x1, x2, x3, x4;
1683  if ( a1x <= a2x ) {
1684  x1 = a1x;
1685  x2 = a2x;
1686  } else {
1687  x1 = a2x;
1688  x2 = a1x;
1689  }
1690  if ( b1x <= b2x ) {
1691  x3 = b1x;
1692  x4 = b2x;
1693  } else {
1694  x3 = b2x;
1695  x4 = b1x;
1696  }
1697  double y1, y2, y3, y4;
1698  if ( a1y <= a2y ) {
1699  y1 = a1y;
1700  y2 = a2y;
1701  } else {
1702  y1 = a2y;
1703  y2 = a1y;
1704  }
1705  if ( b1y <= b2y ) {
1706  y3 = b1y;
1707  y4 = b2y;
1708  } else {
1709  y3 = b2y;
1710  y4 = b1y;
1711  }
1712 
1713  // test
1714  double EPS = 1e-10;
1715  if ( ( x1 <= x + EPS ) && ( x <= x2 + EPS ) && ( x3 <= x + EPS ) && ( x <= x4 + EPS ) && ( y1 <= y + EPS ) && ( y <= y2 + EPS ) && ( y3 <= y + EPS ) && ( y <= y4 + EPS ) ) {
1716  intersectCoords.resize(2);
1717  intersectCoords.at(1) = x;
1718  intersectCoords.at(2) = y;
1719  return true;
1720  } else {
1721  intersectCoords.clear();
1722  return false;
1723  }
1724 }
1725 
1726 
1727 
1728 
1729 void
1731  // transform stiff. tensor 9x9 to stiff. matrix 6x6 (universal for 3D and 2D)
1732  IntArray indicesI, indicesJ, indicesK, indicesL;
1733  indicesI = {
1734  1, 2, 3, 2, 1, 1
1735  };
1736  indicesJ = {
1737  1, 2, 3, 3, 3, 2
1738  };
1739  indicesK = {
1740  1, 2, 3, 2, 1, 1
1741  };
1742  indicesL = {
1743  1, 2, 3, 3, 3, 2
1744  };
1745  int ii, jj, i, j, k, l;
1746 
1747  for ( ii = 1; ii <= 6; ii++ ) {
1748  for ( jj = 1; jj <= 6; jj++ ) {
1749  i = indicesI.at(ii);
1750  j = indicesJ.at(ii);
1751 
1752  k = indicesK.at(jj);
1753  l = indicesL.at(jj);
1754 
1755  matrix->at(ii, jj) = tensor->at(3 * ( i - 1 ) + k, 3 * ( j - 1 ) + l);
1756  }
1757  }
1758 }
1759 } // end namespace oofem
bool contains(int value) const
Definition: intarray.h:283
CrossSection * giveCrossSection()
Definition: element.C:495
void applyApproach3(Domain *d, int homMtrxType)
void computeStiffnessTensorOf1Link(FloatMatrix &D1, double &S0, Element *e, Domain *d)
Material * createMaterial(const char *name, int num, Domain *domain)
Creates new instance of material corresponding to given keyword.
Definition: classfactory.C:199
void setField(int item, InputFieldType id)
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: crosssection.C:67
IntArray interpolationElementNumbers
std::vector< IntArray > interpolationMeshNodes
void setCrossSection(int i, CrossSection *obj)
Sets i-th component. The component will be further managed and maintained by domain object...
Definition: domain.C:456
bool computeIntersectionsOfLinkWith2DTringleElements(IntArray &intersected, FloatArray &lengths, Domain *d, Element *e, qcNode *qn1, qcNode *qn2)
bool computeIntersectionsOfLinkWith3DTetrahedraElements(IntArray &intersected, FloatArray &lengths, Domain *d, Element *e, qcNode *qn1, qcNode *qn2)
void erase(int pos)
Erase the element at given position (1-based index) Receiver will shrink accordingly, the values at positions (pos+1,...,size) will be moved to positions (pos,...,size-1)
Definition: intarray.C:163
Class and object Domain.
Definition: domain.h:115
void transformStiffnessTensorToMatrix(FloatMatrix *matrix, FloatMatrix *tensor)
void applyApproach1(Domain *d)
int giveNumberOfDofManagers() const
Returns number of dof managers in domain.
Definition: domain.h:432
void zero()
Sets all component to zero.
Definition: intarray.C:52
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
int giveGlobalNumber() const
Definition: element.h:1059
void setActivatedNodeList(IntArray nodeList, Domain *d)
int intersectionTestSegmentTriangle2D(FloatArray &intersectCoordsX, FloatArray &intersectCoordsY, FloatArray *A, FloatArray *B, FloatArray *C, FloatArray *U1, FloatArray *U2)
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
virtual FloatArray * giveCoordinates()
Definition: dofmanager.h:382
EngngModel * giveEngngModel()
Returns engineering model to which receiver is associated.
Definition: domain.C:433
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of cross section property.
Abstract base class for all finite elements.
Definition: element.h:145
virtual int giveMasterElementNumber()
Definition: qcnode.h:99
void createInterpolationElements(Domain *d)
CrossSection * createCrossSection(const char *name, int num, Domain *domain)
Creates new instance of cross section corresponding to given keyword.
Definition: classfactory.C:189
#define _IFT_SimpleCrossSection_width
#define S(p)
Definition: mdm.C:481
virtual double giveQcPlasticParamneter()=0
void resizeElements(int _newSize)
Resizes the internal data structure to accommodate space for _newSize elements.
Definition: domain.C:445
int giveNumberOfElements() const
Returns number of elements in domain.
Definition: domain.h:434
int giveNumberOfSpatialDimensions()
Returns number of spatial dimensions.
Definition: domain.C:1067
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
void computeIntersectionsOfLinkWithInterpElements(IntArray &intersected, FloatArray &lengths, Domain *d, Element *e, qcNode *qn1, qcNode *qn2)
virtual int giveNumberOfDofManagers() const
Definition: element.h:656
virtual double giveQcElasticParamneter()=0
void setupInterpolationMesh(Domain *d, int generateInterpolationElements, int interpolationElementsMaterialNumber, std::vector< IntArray > *newMeshNodes)
bool stiffnessAssignment(std::vector< FloatMatrix * > &individualStiffnessTensors, FloatArray &individialS0, Domain *d, Element *e, qcNode *qn1, qcNode *qn2)
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
Base abstract class representing cross section in finite element mesh.
Definition: crosssection.h:107
virtual int giveNumberOfNodes() const
Returns number of nodes of receiver.
Definition: element.h:662
This class implements linear static engineering problem.
void setMaterial(int i, Material *obj)
Sets i-th component. The component will be further managed and maintained by domain object...
Definition: domain.C:457
#define _IFT_IsotropicLinearElasticMaterial_talpha
#define _IFT_IsotropicLinearElasticMaterial_n
const IntArray & giveDofManArray() const
Definition: element.h:592
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: element.C:638
Element * giveElement(int n)
Service for accessing particular domain fe element.
Definition: domain.C:160
bool intersectionTestSegmentTrianglePlucker3D(FloatArray &intersectCoords, FloatArray *A, FloatArray *B, FloatArray *C, FloatArray *X1, FloatArray *X2)
void homogenizationOfStiffMatrix(double &homogenizedE, double &homogenizedNu, FloatMatrix *Diso)
#define OOFEM_ERROR(...)
Definition: error.h:61
void clear()
Clears the array (zero size).
Definition: intarray.h:177
void setActivatedElementList(IntArray elemList)
#define _IFT_Material_density
Definition: material.h:51
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
#define _IFT_Element_nodes
Definition: element.h:65
Class implementing hanging node connected to other nodes (masters) using interpolation.
Definition: qcnode.h:64
int giveNumberOfCrossSectionModels() const
Returns number of cross section models in domain.
Definition: domain.h:438
int giveNumberOfMaterialModels() const
Returns number of material models in domain.
Definition: domain.h:436
virtual void setCrossSection(int csIndx)
Sets the cross section model of receiver.
Definition: element.h:653
#define _IFT_IsotropicLinearElasticMaterial_e
void setNoDimensions(Domain *d)
Abstract base class for all material models.
Definition: material.h:95
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
const IntArray & giveElementsWithMaterialNum(int iMaterialNum) const
Returns array with indices of elements that have a given material number.
Definition: domain.C:209
Material interface for gradient material models.
int giveIndexMaxElem(void)
Returns index (between 1 and Size) of maximum element in the array.
Definition: floatarray.C:447
void createGlobalStiffnesMatrix(FloatMatrix *Diso, double &S0, Domain *d, int homMtrxType, double volumeOfInterpolationMesh)
void applyApproach2(Domain *d, int homMtrxType, double volumeOfInterpolationMesh)
Class representing vector of real numbers.
Definition: floatarray.h:82
void resizeMaterials(int _newSize)
Resizes the internal data structure to accommodate space for _newSize materials.
Definition: domain.C:447
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
Class implementing "simple" cross section model in finite element problem.
#define _IFT_AnisotropicLinearElasticMaterial_talpha
void initializeConnectivityTableForInterpolationElements(Domain *d)
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
void addCrosssectionToInterpolationElements(Domain *d)
void add(const FloatMatrix &a)
Adds matrix to the receiver.
Definition: floatmatrix.C:1023
void followedBy(const IntArray &b, int allocChunk=0)
Appends array b at the end of receiver.
Definition: intarray.C:145
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual int giveQcNodeType()
Definition: node.h:183
void setMaterial(int matIndx)
Sets the material of receiver.
Definition: element.h:647
void setGlobalNumber(int num)
Sets receiver globally unique number.
Definition: element.h:1064
void setElement(int i, Element *obj)
Sets i-th component. The component will be further managed and maintained by domain object...
Definition: domain.C:455
virtual Interface * giveInterface(InterfaceType t)
Interface requesting service.
Definition: femcmpnn.h:179
Class representing the a dynamic Input Record.
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
ClassFactory & classFactory
Definition: classfactory.C:59
int intersectionTestSegmentTetrahedra3D(FloatArray &intersectCoordsX, FloatArray &intersectCoordsY, FloatArray &intersectCoordsZ, FloatArray *A, FloatArray *B, FloatArray *C, FloatArray *D, FloatArray *X1, FloatArray *X2)
virtual FloatArray * giveCoordinates()
Definition: node.h:114
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
void push_back(const double &iVal)
Add one element.
Definition: floatarray.h:118
std::vector< IntArray * > connectivityTable
void resizeCrossSectionModels(int _newSize)
Resizes the internal data structure to accommodate space for _newSize cross section models...
Definition: domain.C:446
#define _IFT_SimpleCrossSection_thick
bool intersectionTestSegmentSegment2D(FloatArray &intersectCoords, FloatArray *A1, FloatArray *A2, FloatArray *B1, FloatArray *B2)
int giveSize() const
Definition: intarray.h:203
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
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
#define _IFT_AnisotropicLinearElasticMaterial_stiff
int giveNumber() const
Definition: femcmpnn.h:107
Element * createElement(const char *name, int num, Domain *domain)
Creates new instance of element corresponding to given keyword.
Definition: classfactory.C:159
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: material.C:89
IntArray interpolationElementIndices
#define OOFEM_WARNING(...)
Definition: error.h:62
int findFirstIndexOf(int value) const
Finds index of first occurrence of given value in array.
Definition: intarray.C:331
virtual Material * giveMaterial()
Definition: element.C:484
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631

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:31 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011