OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
unstructuredgridfield.h
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 
35 #ifndef unstructuredgridfield_h
36 #define unstructuredgridfield_h
37 
38 #include "field.h"
39 #include "floatarray.h"
40 #include "intarray.h"
41 #include "elementgeometrytype.h"
42 #include "feinterpol.h"
43 #include "fei2dlinelin.h"
44 #include "fei2dlinequad.h"
45 #include "fei2dtrlin.h"
46 #include "fei2dtrquad.h"
47 #include "fei2dquadlin.h"
48 #include "fei2dquadquad.h"
49 #include "fei3dtetlin.h"
50 #include "fei3dhexalin.h"
51 
52 #include "pfem/octreelocalizert.h"
53 #include "error.h"
54 
55 namespace oofem {
56 
63 class OOFEM_EXPORT UnstructuredGridField: public Field
64 {
65 protected:
66 
67 
68 
69  class Vertex {
71  public:
72  Vertex (): coords() {}
73  Vertex (FloatArray& c) {coords = c;}
74  Vertex& operator= (const Vertex& v) {
75  coords = v.coords;
76  return *this;
77  }
78  const FloatArray* getCoordinates() const {return &coords;}
79  };
80 
81 
82 
83 
84  /* Auxiliary class to represent interpolation cell, basically holding vertices defining its shape and interpolation type */
85  class Cell {
89 
90  // array of interpolation instances used for supported ElementGeometryTypes
91  static FEI2dLineLin i1;
92  static FEI2dLineQuad i2;
93  static FEI2dTrLin i3;
94  static FEI2dTrQuad i4;
95  static FEI2dQuadLin i5;
96  static FEI2dQuadQuad i6;
97  static FEI3dTetLin i7;
98  static FEI3dHexaLin i8;
99 
100  static FEInterpolation* interpTable[] ;
101 
102 
107  {
108  protected:
109  const Cell *cell;
110  public:
112  this->cell = c;
113  }
115  int giveNumberOfVertices() const {
116  return this->cell->giveNumberOfVertices();
117  }
118  inline const FloatArray *giveVertexCoordinates(int i) const
119  {
120  return ((cell->getVertex(i-1))->getCoordinates());
121  }
122  };
123 
124  public:
125  Cell () {
126  itype = EGT_unknown;
127  mesh=NULL;
128  }
130  itype = t;
131  vertices = v;
132  mesh = m;
133  }
134 
135  Cell& operator= (const Cell& c) {
136  itype = c.itype;
137  vertices = c.vertices;
138  mesh = c.mesh;
139  return *this;
140  }
141 
142  int giveNumberOfVertices() const {return vertices.giveSize();}
143  bool containsPoint(const FloatArray &coords) const {
144  FloatArray tmp;
145  BoundingBox b;
146  this->giveBoundingBox(b);
147  if (b.contains(coords)) {
149  return (this->getInterpolation()->global2local (tmp, coords, cgw));
150  }
151  return false;
152  }
153 
154  void giveBoundingBox(BoundingBox& bb) const {
155  double size;
156  FloatArray bb0, bb1;
157  bb1 = bb0 = *(this->getVertex(0)->getCoordinates());
158 
159  for ( int i = 1; i < this->giveNumberOfVertices(); ++i ) {
160  const FloatArray *coordinates = this->getVertex(i)->getCoordinates();
161  bb0.beMinOf(bb0, * coordinates);
162  bb1.beMaxOf(bb1, * coordinates);
163  }
164  bb1.subtract(bb0);
165  int nsd = bb1.giveSize();
166  IntArray mask(3);
167  size = bb1.at(1);
168  for (int i =1; i< nsd; i++) size=max(size, bb1(i));
169  for (int i=0; i<nsd; i++) mask(i)=1;
170  bb.init (bb0, size, mask);
171  }
172 
173  double giveClosestPoint(FloatArray &lcoords, FloatArray &closest, const FloatArray &gcoords) {
174  FEInterpolation *interp = this->getInterpolation();
175 
176  if ( !interp->global2local( lcoords, gcoords, FEICellGeometryWrapper(this) ) ) { // Outside element
177  interp->local2global( closest, lcoords, FEICellGeometryWrapper(this) );
178  return closest.distance(gcoords);
179  } else {
180  closest = gcoords;
181  return 0.0;
182  }
183  }
184 
186  if (this->itype == EGT_line_1) return interpTable[0];
187  else if (this->itype == EGT_line_2) return interpTable[1];
188  else if (this->itype == EGT_triangle_1) return interpTable[2];
189  else if (this->itype == EGT_triangle_2) return interpTable[3];
190  else if (this->itype == EGT_quad_1) return interpTable[4];
191  else if (this->itype == EGT_quad_2) return interpTable[5];
192  else if (this->itype == EGT_tetra_1) return interpTable[6];
193  else if (this->itype == EGT_hexa_1) return interpTable[7];
194  else {
195  OOFEM_LOG_ERROR ("UnstructuredGridField.Cell:: Unsupported cell type");
196  return NULL;
197  }
198 
199  }
200  const Vertex* getVertex (int i) const {return mesh->getVertex(vertices(i));}
201  int interpolate (FloatArray& answer, const FloatArray& pos, FloatArray** vertexVals) {
202  int i,j, size = vertexVals[0]->giveSize();
203  FloatArray N, lcoords;
204 
205  FEInterpolation* it = this->getInterpolation();
206  FEICellGeometryWrapper cw (this);
207  it->global2local(lcoords, pos, cw);
208  it->evalN (N, lcoords, cw);
209 
210  answer.resize(size);
211  answer.zero();
212  for (i=0; i<giveNumberOfVertices(); i++) {
213  for (j=0; j<size; j++) {
214  answer(j)+=N(i)*vertexVals[i]->at(j+1);
215  }
216  }
217  return 1;
218  }
219  int getVertexNum (int i) {
220  return this->vertices(i);
221  }
222  };
223 
224 
226  private:
227  public:
229  virtual bool evaluate (Cell& member, OctantRecT<Cell>* cell) {
230  BoundingBox b;
231  member.giveBoundingBox(b);
234  return true;
235  } else {
236  return false;
237  }
238  }
239 
241  std::list<LocalInsertionData<Cell>> *giveInsertionList(Cell& m) {return NULL;}
242 
243  };
244 
246  {
247  protected:
249  std :: list <Cell > cells;
250  public:
256  CellContainingPointFunctor(const FloatArray &pos) : position (pos) {}
258 
264  bool evaluate(Cell& c)
265  {
266  if (c.containsPoint(position)) {
267  cells.push_back(c);
268  return true;
269  } else {
270  return false;
271  }
272  }
273 
279  {
280  answer = position;
281  }
282 
287  void giveResult(std :: list <Cell > &answer)
288  { answer=cells;}
289 
290  bool isBBXStage1Defined(BoundingBox &BBXStage1) { return true; }
291  bool isBBXStage2Defined(BoundingBox &BBXStage2) { return false; }
292 
293  };
294 
295 
296  std::vector<Cell> cellList;
297  std::vector<Vertex> vertexList;
298  std::vector<FloatArray> valueList;
305  long int octreeTimeStamp;
307  long int timeStamp;
310  public:
314  UnstructuredGridField(int nvert, int ncells, double octreeOriginShift=0.0 ) : Field(FieldType::FT_Unknown), spatialLocalizer()
315  { this->timeStamp = this->octreeTimeStamp = 0;
316  this->vertexList.resize(nvert);
317  this->cellList.resize(ncells);
318  this->valueList.resize(nvert);
319  this->octreeOriginShift = octreeOriginShift;
320  }
322 
323  void setVertexValue(int num, const FloatArray& vv) {
324  valueList[num] = vv;
325  }
326 
327  void addVertex (int num, FloatArray& coords) {
328  vertexList[num] = Vertex(coords);
329  this->timeStamp++;
330  }
331 
332  Vertex* getVertex (int num) {return &this->vertexList[num];}
333 
334  void addCell (int num, Element_Geometry_Type type, IntArray& vertices) {
335  cellList[num]=Cell(type, vertices, this);
336  this->timeStamp ++;
337  }
338 
339  int evaluateAt(FloatArray &answer, const FloatArray &coords,
340  ValueModeType mode, TimeStep *tStep) override {
341  std::list<Cell> elist;
342  if ((mode == VM_Total) || (mode == VM_TotalIntrinsic)) {
343  if(this->cellList.size()>0){
344  this->initOctree();
345  CellContainingPointFunctor f(coords);
346  this->spatialLocalizer.giveDataOnFilter(elist, f);
347  if (elist.size()) {
348  Cell &c = elist.front(); // take first
349  // colect vertex values
350  int size = c.giveNumberOfVertices();
351  FloatArray** vertexValues = new FloatArray*[size];
352  for (int i=0; i<size; i++) {
353  vertexValues[i]=&(this->valueList[c.getVertexNum(i)]);
354  }
355  c.interpolate (answer, coords, vertexValues);
356  return 0;
357  } else {
358  return 1;
359  }
360  } else {
361  if(!this->vertexList.size()) {
362  return 1;
363  }
364  // alternative method - we use the value of the closest point
365  double minDist=0.,dist=0.;
366 
367  int idOfClosestPoint=-1;
368  for(int i=0;i<(int) this->vertexList.size();i++){
369  const FloatArray *pcoords=this->vertexList[i].getCoordinates();
370  dist=sqrt( pow(coords[0]-pcoords->at(1),2)+pow(coords[1]-pcoords->at(2),2)+pow(coords[2]-pcoords->at(3),2) );
371  if((dist<minDist) || (!i)){
372  minDist=dist;
373  idOfClosestPoint=i;
374  }
375  }
376  answer=this->valueList[idOfClosestPoint];
377  //printf("closest point is %d\n",idOfClosestPoint);
378  return 0;
379  }
380  } else {
381  OOFEM_ERROR("Unsupported ValueModeType");
382  }
383  }
384 
388  int evaluateAt(FloatArray &answer, DofManager *dman,
389  ValueModeType mode, TimeStep *tStep) override {
390  FloatArray* coords = dman->giveCoordinates();
391  return this->evaluateAt (answer, *coords, mode, tStep);
392  }
393 
402  contextIOResultType saveContext(DataStream &stream, ContextMode mode) override { return CIO_OK; }
411  contextIOResultType restoreContext(DataStream &stream, ContextMode mode) override { return CIO_OK; }
412 
414  const char *giveClassName() const override { return "UnstructuredGridField"; }
415 
416  protected:
417  void initOctree() {
418  if (this->timeStamp != this->octreeTimeStamp) {
419  // rebuild octree
420  this->spatialLocalizer.clear();
421  // get octree bbox
422  std::vector<Vertex>::iterator it = vertexList.begin();
423  FloatArray cmax, cmin;
424  cmax= cmin = *((*it).getCoordinates());
425  int nsd = cmin.giveSize();
426  ++it;
427  for (; it != vertexList.end(); ++it) {
428  const FloatArray *vc = (*it).getCoordinates();
429  for (int j=0; j<nsd; j++) {
430  cmax(j)=max(cmax(j), (*vc)(j));
431  cmin(j)=min(cmin(j), (*vc)(j));
432  }
433  }
434  // introduce origin shift (if set up)
435  for (int j=0; j<nsd; j++) {
436  cmin(j) -= octreeOriginShift;
437  }
438  double size = 0.0;
439  IntArray mask (3);
440  for (int j=0; j<nsd; j++) {
441  size = max(size, cmax(j)-cmin(j));
442  mask(j)=1;
443  }
444  BoundingBox bb;
445  bb.init(cmin, size, mask);
446  this->spatialLocalizer.init (bb);
447  std::vector<Cell>::iterator cit;
449  //int cellcount = 1;
450  for (cit=this->cellList.begin(); cit != this->cellList.end(); ++cit) {
451  /*
452  BoundingBox b; (*cit).giveBoundingBox(b);
453  FloatArray o;
454  b.giveOrigin(o);
455  printf ("Inserting cell %d, Origin, Size:", cellcount++);
456  o.printYourself();
457  printf("%lf\n", b.giveSize());
458  */
459  this->spatialLocalizer.insertMemberIntoOctree(*cit, cf);
460  }
461 
462  // update timestamp
463  this->octreeTimeStamp = this->timeStamp;
464  }
465  }
466 };
467 } // end namespace oofem
468 #endif // unstructuredgridfield_h
Class representing implementation of linear hexahedra interpolation class.
Definition: fei3dhexalin.h:44
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
Templated octree cell containing data of T type.
bool isBBXStage1Defined(BoundingBox &BBXStage1)
Stage1 means, we are looking for objects in a distance given by some boundingBox (e.g.
Squared bounding box for templated octree localizer.
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)=0
Evaluates the array of interpolation functions (shape functions) at given point.
#define OOFEM_LOG_ERROR(...)
Definition: logger.h:122
FieldType
Physical type of field.
Definition: field.h:60
Element_Geometry_Type
Enumerative type used to classify element geometry Possible values are: EGT_point - point in space EG...
const Vertex * getVertex(int i) const
int interpolate(FloatArray &answer, const FloatArray &pos, FloatArray **vertexVals)
std::list< LocalInsertionData< Cell > > * giveInsertionList(Cell &m)
const FloatArray * getCoordinates() const
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
int max(int i, int j)
Returns bigger value form two given decimals.
Definition: mathfem.h:71
ValueModeType
Type representing the mode of UnknownType or CharType, or similar types.
Definition: valuemodetype.h:78
Class representing a general abstraction for cell geometry.
Definition: feinterpol.h:62
Class representing implementation of linear tetrahedra interpolation class.
Definition: fei3dtetlin.h:44
contextIOResultType restoreContext(DataStream &stream, ContextMode mode) override
Restores the receiver state previously written in stream.
bool evaluate(Cell &c)
Evaluates a triangle upon its circumscribed cricle.
virtual FloatArray * giveCoordinates()
Definition: dofmanager.h:382
Cell(Element_Geometry_Type t, IntArray &v, UnstructuredGridField *m)
long int octreeTimeStamp
octree build time stamp
void init(FloatArray &origin, double size, IntArray &mask)
Sets all BBOx parameters in ince.
Base class for dof managers.
Definition: dofmanager.h:113
std::vector< FloatArray > valueList
void giveResult(std::list< Cell > &answer)
Gives the triangles containing the node.
void giveDataOnFilter(std::list< T > &answer, SL_Evaluation_Functor< T > &filter)
Evalutes the search accoring used functor a fills the list with results - NOT IN USE.
Class implementing an array of integers.
Definition: intarray.h:61
Wrapper around element definition to provide FEICellGeometry interface.
Class representing a 2d line quadratic interpolation.
Definition: fei2dlinequad.h:46
void beMaxOf(const FloatArray &a, const FloatArray &b)
Sets receiver to maximum of a or b&#39;s respective elements.
Definition: floatarray.C:288
Abstract class representing field.
Definition: field.h:80
double distance(const FloatArray &x) const
Computes the distance between position represented by receiver and position given as parameter...
Definition: floatarray.C:489
Class representing a general abstraction for finite element interpolation class.
Definition: feinterpol.h:132
UnstructuredGridField(int nvert, int ncells, double octreeOriginShift=0.0)
Constructor.
Class representing a 2d triangular linear interpolation based on area coordinates.
Definition: fei2dtrlin.h:44
FEInterpolation * getInterpolation() const
void beMinOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be minimum of a or b&#39;s respective elements.
Definition: floatarray.C:315
int evaluateAt(FloatArray &answer, const FloatArray &coords, ValueModeType mode, TimeStep *tStep) override
Evaluates the field at given point.
virtual bool evaluate(Cell &member, OctantRecT< Cell > *cell)
#define OOFEM_ERROR(...)
Definition: error.h:61
void addVertex(int num, FloatArray &coords)
Help class for storing pointer to octant cell and position of the member in the data list...
virtual int global2local(FloatArray &answer, const FloatArray &gcoords, const FEICellGeometry &cellgeo)=0
Evaluates local coordinates from given global ones.
#define N(p, q)
Definition: mdm.C:367
Field defined by values fefined on unstructured grid.
Class representing a 2d quadrilateral with quadratic interpolation based on isoparametric coordinates...
Definition: fei2dquadquad.h:52
void addCell(int num, Element_Geometry_Type type, IntArray &vertices)
void setVertexValue(int num, const FloatArray &vv)
Class representing a 2d line with linear interpolation.
Definition: fei2dlinelin.h:46
void giveBoundingBox(BoundingBox &bb) const
Class representing vector of real numbers.
Definition: floatarray.h:82
bool isBBXStage2Defined(BoundingBox &BBXStage2)
Stage2BBX is given by results of a prior search.
contextIOResultType saveContext(DataStream &stream, ContextMode mode) override
Stores receiver state to output stream.
Functor base class responsible for insertion of members into the octree cell.
bool contains(const FloatArray &coords) const
OctantRec::BoundingBoxStatus testBoundingBox(BoundingBox &testedBBX)
Tests the position of a bounding box in relation to the octant cell.
bool containsPoint(const FloatArray &coords) const
const char * giveClassName() const override
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
int insertMemberIntoOctree(T &memberID, SL_Insertion_Functor< T > &functor)
Inserts member into the octree using functor for the evaluation.
double octreeOriginShift
octree origin shift
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
Templated octree spatial localizer.
Functor base class for evaluating search tasks on the octree according given condition.
CellContainingPointFunctor(const FloatArray &pos)
Constructor.
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
int evaluateAt(FloatArray &answer, DofManager *dman, ValueModeType mode, TimeStep *tStep) override
Implementaton of Field::evaluateAt for DofManager.
int init(BoundingBox &BBX, int initialDivision=0)
Initilizes the octree structure.
long int timeStamp
receiver timestamp
Second order triangular interpolation in 2D (6 nodes).
Definition: fei2dtrquad.h:44
Class representing a 2d isoparametric linear interpolation based on natural coordinates for quadrilat...
Definition: fei2dquadlin.h:45
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.
double giveClosestPoint(FloatArray &lcoords, FloatArray &closest, const FloatArray &gcoords)
void giveStartingPosition(FloatArray &answer)
Gives the starting position of the search.
OctreeSpatialLocalizerT< Cell > spatialLocalizer
Spatial Localizer.
Class representing solution step.
Definition: timestep.h:80
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)=0
Evaluates global coordinates from given local ones.
void registerInsertion(Cell &member, LocalInsertionData< Cell > lidata)
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:32 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011