geometry.cpp 44.2 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#include "inmost.h"
#if defined(USE_MESH)


namespace INMOST
{
	
	__INLINE static void vec_diff(const Storage::real * vecin1, const Storage::real * vecin2, Storage::real * vecout, unsigned int size)
	{
		for(unsigned int i = 0; i < size; i++)
			vecout[i] = vecin1[i] - vecin2[i];
	}
	
	__INLINE static void vec_cross_product(const Storage::real * vecin1, const Storage::real * vecin2, Storage::real * vecout)
	{
		Storage::real temp[3];
		temp[0] = vecin1[1]*vecin2[2] - vecin1[2]*vecin2[1];
		temp[1] = vecin1[2]*vecin2[0] - vecin1[0]*vecin2[2];
		temp[2] = vecin1[0]*vecin2[1] - vecin1[1]*vecin2[0];
		vecout[0] = temp[0];
		vecout[1] = temp[1];
		vecout[2] = temp[2];
	}
	
	__INLINE static Storage::real vec_dot_product(const Storage::real * vecin1,const Storage::real * vecin2, unsigned int size)
	{
		Storage::real ret = 0;
		for(unsigned int i = 0; i < size; i++)
			ret += vecin1[i]*vecin2[i];
		return ret;
	}
	
	__INLINE static Storage::real vec_len2(const Storage::real * vecin, unsigned int size)
	{
		return vec_dot_product(vecin,vecin,size);
	}
	
	__INLINE static Storage::real vec_len(const Storage::real * vecin, unsigned int size)
	{
		return sqrt(vec_len2(vecin,size));
	}
	
	__INLINE static Storage::real det3d(Storage::real a, Storage::real b, Storage::real c,
	                         Storage::real d, Storage::real e, Storage::real f,
	                         Storage::real g, Storage::real h, Storage::real i ) 
	{
		return a*e*i - c*e*g + b*f*g - a*f*h + c*d*h - b*d*i;
	}
	
	__INLINE static Storage::real det3v(const Storage::real * x,const Storage::real * y,const Storage::real * z) 
	{
		return det3d(x[0], x[1], x[2],  y[0], y[1], y[2],  z[0], z[1], z[2]);
	}
	
	__INLINE static Storage::real det4v(const Storage::real * w, const Storage::real * x, const Storage::real * y, const Storage::real * z) 
	{
		return det3d(x[0]-w[0], x[1]-w[1], x[2]-w[2],  y[0]-w[0], y[1]-w[1], y[2]-w[2],  z[0]-w[0], z[1]-w[1], z[2]-w[2]);
	}	
	
	
	__INLINE static Storage::real vec_normalize(Storage::real * vecin, unsigned int size)
	{
		Storage::real len = 0;
		for(unsigned int i = 0; i < size; i++)
			len += vecin[i]*vecin[i];
		len = sqrt(len);
		for(unsigned int i = 0; i < size; i++)
			vecin[i] /= len;
		return len;
	}
	
	
Kirill Terekhov's avatar
Kirill Terekhov committed
73
	ElementArray<Cell> Cell::NeighbouringCells() const
Kirill Terekhov's avatar
Kirill Terekhov committed
74
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
75
76
77
		ElementArray<Cell> ret(GetMeshLink());
		ElementArray<Face> faces = getFaces();
		for(ElementArray<Face>::iterator f = faces.begin(); f != faces.end(); f++)
Kirill Terekhov's avatar
Kirill Terekhov committed
78
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
79
80
			Cell c = Neighbour(f->self());
			if( c.isValid() ) ret.push_back(c);
Kirill Terekhov's avatar
Kirill Terekhov committed
81
82
83
84
85
		}
		return ret;
	}
	
	
Kirill Terekhov's avatar
Kirill Terekhov committed
86
	Cell Cell::Neighbour(Face f) const
Kirill Terekhov's avatar
Kirill Terekhov committed
87
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
88
89
		Cell b = f->BackCell();
		if( b == self() )
Kirill Terekhov's avatar
Kirill Terekhov committed
90
91
92
93
			return f->FrontCell();
		return b;
	}
	
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
94
	bool Cell::Inside(const Storage::real * point) const//check for 2d case
Kirill Terekhov's avatar
Kirill Terekhov committed
95
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
96
		Mesh * m = GetMeshLink();
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
		integer dim = GetElementDimension();
		if( dim == 3 )
		{
			assert(m->GetDimensions() == 3);
			integer vp = 0;
			integer vm = 0;
			integer vz = 0;
			real eps = m->GetEpsilon();
			real c,d;
			adj_type & elements = m->LowConn(GetHandle());
			for(adj_type::size_type f = 0; f < elements.size(); f++)
			{
				Face ef = Face(m,elements[f]);
				d = 0.0;
				ElementArray<Node> nodes = ef->getNodes();
				for (ElementArray<Node>::size_type i=1; i<nodes.size()-1; i++) 
					d += c = det4v(point, nodes[0].Coords().data(), nodes[i].Coords().data(), nodes[i+1].Coords().data());
				if(ef->FaceOrientedOutside(self()) == 0)  
					c = -1.0;  
				else 
					c = 1.0;
				if(c*d > eps)  
					vp++;  
				else if(c*d < eps)  
					vm++;  
				else  
					vz++;
			}
			if(vp*vm > 0) return false;
			return false;
		}
		else
Kirill Terekhov's avatar
Kirill Terekhov committed
129
		{
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
130
131
132
133
134
135
136
137
138
			int mdim = m->GetDimensions();
			assert(mdim <= 3);
			Storage::real data[9][3];
			if( mdim < 3 )
			{
				memset(data,0,sizeof(Storage::real)*9*3);
			}
			Centroid(data[0]);
			ElementArray<Node> nodes = getNodes();
Kirill Terekhov's avatar
Kirill Terekhov committed
139
			for(int i = 0; i < static_cast<int>(nodes.size()); i++)
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
			{
				int j = (i+1)%nodes.size();
				nodes[i].Centroid(data[1]);
				nodes[j].Centroid(data[2]);
				vec_diff(point,data[0],data[3],mdim);
				vec_diff(point,data[1],data[4],mdim);
				vec_diff(point,data[2],data[5],mdim);
				vec_cross_product(data[3],data[4],data[6]);
				vec_cross_product(data[4],data[5],data[7]);
				vec_cross_product(data[5],data[3],data[8]);

				if( vec_dot_product(data[6],data[7],mdim) >= 0 &&
				    vec_dot_product(data[7],data[8],mdim) >= 0 &&
				    vec_dot_product(data[6],data[8],mdim) >= 0 )
					return true; //inside one of the triangles
			}
			return false;
Kirill Terekhov's avatar
Kirill Terekhov committed
157
158
159
160
		}
	}
	
	
Kirill Terekhov's avatar
Kirill Terekhov committed
161
	void Face::UnitNormal(real * nrm) const
Kirill Terekhov's avatar
Kirill Terekhov committed
162
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
163
164
165
166
167
168
169
		Mesh * m = GetMeshLink();
		m->GetGeometricData(GetHandle(),NORMAL,nrm); 
		integer dim = m->GetDimensions();
		real    l   = sqrt(vec_dot_product(nrm,nrm,dim)); 
		if(fabs(l) > m->GetEpsilon()) 
		{
			for(integer i = 0; i < dim; i++) 
Kirill Terekhov's avatar
Kirill Terekhov committed
170
				nrm[i] /= l; 
Kirill Terekhov's avatar
Kirill Terekhov committed
171
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
172
173
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
174
	void Face::OrientedNormal(Cell c, Storage::real * nrm) const
Kirill Terekhov's avatar
Kirill Terekhov committed
175
176
	{
		Normal(nrm); 
Kirill Terekhov's avatar
Kirill Terekhov committed
177
178
179
180
		if( !FaceOrientedOutside(c) )
		{
			integer dim = GetMeshLink()->GetDimensions();
			for(integer i = 0; i < dim; i++) 
Kirill Terekhov's avatar
Kirill Terekhov committed
181
				nrm[i] = -nrm[i];
Kirill Terekhov's avatar
Kirill Terekhov committed
182
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
183
184
185
	}
	
	
Kirill Terekhov's avatar
Kirill Terekhov committed
186
	void Face::OrientedUnitNormal(Cell c, Storage::real * nrm) const
Kirill Terekhov's avatar
Kirill Terekhov committed
187
188
189
	{
		UnitNormal(nrm); 
		if( !FaceOrientedOutside(c) ) 
Kirill Terekhov's avatar
Kirill Terekhov committed
190
191
192
		{
			integer dim = GetMeshLink()->GetDimensions();
			for(integer i = 0; i < dim; i++) 
Kirill Terekhov's avatar
Kirill Terekhov committed
193
				nrm[i] = -nrm[i];
Kirill Terekhov's avatar
Kirill Terekhov committed
194
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
195
196
197
198
	}
	
	
	
Kirill Terekhov's avatar
Kirill Terekhov committed
199
	bool Mesh::TestClosure(const HandleType * elements, integer num) const
Kirill Terekhov's avatar
Kirill Terekhov committed
200
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
201
202
203
		integer i;
		tiny_map<HandleType,int,64> e_visit;
		tiny_map<HandleType,int,64>::iterator it;
Kirill Terekhov's avatar
Kirill Terekhov committed
204
205
206
207
		if( !HideMarker() )
		{
			for(i = 0; i < num; i++)
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
208
209
210
				Element::adj_type const & lc = LowConn(elements[i]);
				for(Element::adj_type::size_type jt = 0; jt < lc.size(); jt++)
					e_visit[lc[jt]]++;
Kirill Terekhov's avatar
Kirill Terekhov committed
211
212
213
214
			}
		}
		else
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
215
			for(i = 0; i < num; i++) if( !GetMarker(elements[i],HideMarker()) )
Kirill Terekhov's avatar
Kirill Terekhov committed
216
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
217
218
219
220
221
222
				Element::adj_type const & lc = LowConn(elements[i]);
				for(Element::adj_type::size_type jt = 0; jt < lc.size(); jt++) 
				{
					if( !GetMarker(lc[jt],HideMarker()) ) 
						e_visit[lc[jt]]++;
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
223
224
225
226
227
228
229
			}
		}
		for(it = e_visit.begin(); it != e_visit.end(); it++)
			if( it->second != 2 ) return false;
		return true;
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
230
	Element::GeometricType Mesh::ComputeGeometricType(ElementType etype, const HandleType * lc, INMOST_DATA_ENUM_TYPE s) const
Kirill Terekhov's avatar
Kirill Terekhov committed
231
232
233
234
235
236
237
238
239
240
241
242
243
	{
		Element::GeometricType ret = Element::Unset;
		if( s == 0 && etype != NODE) return ret;
		switch(etype)
		{
			case NODE: ret = Element::Vertex; break;
			case EDGE:
				if( s == 1 )
					ret = Element::Vertex;
				else if( s == 2 )
					ret = Element::Line;
				break;
			case FACE:
Kirill Terekhov's avatar
Kirill Terekhov committed
244
245

				if( Element::GetGeometricDimension(GetGeometricType(lc[0])) == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
				{ 
					ret = Element::Line;
				}
				else
				{
					if( !GetTopologyCheck(NEED_TEST_CLOSURE) || TestClosure(lc,s) )
					{
						if( s == 3 )
							ret = Element::Tri;
						else if( s == 4 )
							ret = Element::Quad;
						else
							ret = Element::Polygon;
					}
					else ret = Element::MultiLine;
				}
				break;
			case CELL:
Kirill Terekhov's avatar
Kirill Terekhov committed
264
				if(  Element::GetGeometricDimension(GetGeometricType(lc[0])) == 1 )
Kirill Terekhov's avatar
Kirill Terekhov committed
265
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
266
					if( !GetTopologyCheck(NEED_TEST_CLOSURE) || TestClosure(lc,s) )
Kirill Terekhov's avatar
Kirill Terekhov committed
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
					{
						if( s == 3 )
							ret = Element::Tri;
						else if( s == 4 )
							ret = Element::Quad;
						else
							ret = Element::Polygon;
					}
					else ret = Element::MultiLine;
				}
				else 
				{
					if( !GetTopologyCheck(NEED_TEST_CLOSURE) ||  TestClosure(lc,s) )
					{
						//test c_faces closure, if no closure, set as MultiPolygon
						INMOST_DATA_ENUM_TYPE quads = 0,tris = 0,i;
						for(i = 0; i < s; i++)
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
285
							if( GetGeometricType(lc[i]) == Element::Tri )
Kirill Terekhov's avatar
Kirill Terekhov committed
286
								tris++;
Kirill Terekhov's avatar
Kirill Terekhov committed
287
							else if( GetGeometricType(lc[i]) == Element::Quad )
Kirill Terekhov's avatar
Kirill Terekhov committed
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
								quads++;
						}
						if( tris == 4 && s == 4 )
							ret = Element::Tet;
						else if( quads == 6 && s == 6 )
							ret = Element::Hex;
						else if( tris == 4 && quads == 1 && s == tris+quads)
							ret = Element::Pyramid;
						else if( quads == 3 && tris == 2 && s == tris+quads)
							ret = Element::Prism;
						else
							ret = Element::Polyhedron;
					}
					else ret = Element::MultiPolygon;
				}
				break;
			case ESET: ret = Element::Set; break;
		}
		return ret;
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
309
	Storage::real Edge::Length() const 
Kirill Terekhov's avatar
Kirill Terekhov committed
310
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
		Storage::real ret; 
		GetMeshLink()->GetGeometricData(GetHandle(),MEASURE,&ret); 
		return ret;
	}

	Storage::real Face::Area() const 
	{
		real ret; 
		GetMeshLink()->GetGeometricData(GetHandle(),MEASURE,&ret); 
		return ret;
	}
	void Face::Normal(real * nrm) const 
	{
		GetMeshLink()->GetGeometricData(GetHandle(),NORMAL,nrm);
	}

	Storage::real Cell::Volume() const 
	{
		real ret; 
		GetMeshLink()->GetGeometricData(GetHandle(),MEASURE,&ret); 
		return ret;
	}

	void Element::ComputeGeometricType() const
	{
		GetMeshLink()->ComputeGeometricType(GetHandle());
Kirill Terekhov's avatar
Kirill Terekhov committed
337
338
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
339
340
341
342
343
344
345
	void Mesh::ComputeGeometricType(HandleType h) 
	{
		SetGeometricType(h,Element::Unset);
		Element::adj_type const & lc = LowConn(h);
		if( !lc.empty() )
			SetGeometricType(h,ComputeGeometricType(GetHandleElementType(h),lc.data(),static_cast<integer>(lc.size())));
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
346
	
Kirill Terekhov's avatar
Kirill Terekhov committed
347
	void Mesh::RecomputeGeometricData(HandleType e)
Kirill Terekhov's avatar
Kirill Terekhov committed
348
349
350
351
352
	{
		//static std::map<Element *, int> numfixes;
		GeometricData d ;
		for(d = CENTROID; d <= NORMAL; d++) // first compute centroids and normals 
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
353
			if( HaveGeometricData(d,GetHandleElementType(e)) ) //compute centroid first
Kirill Terekhov's avatar
Kirill Terekhov committed
354
355
			{
				Tag t = GetGeometricTag(d);
Kirill Terekhov's avatar
Kirill Terekhov committed
356
357
				Storage::real * a = static_cast<Storage::real *>(MGetDenseLink(e,t));
				HideGeometricData(d,GetHandleElementType(e));
Kirill Terekhov's avatar
Kirill Terekhov committed
358
				GetGeometricData(e,d,a);
Kirill Terekhov's avatar
Kirill Terekhov committed
359
				ShowGeometricData(d,GetHandleElementType(e));
Kirill Terekhov's avatar
Kirill Terekhov committed
360
361
362
363
			}
		}


Kirill Terekhov's avatar
Kirill Terekhov committed
364
		if( GetHandleElementType(e) == CELL && HaveGeometricData(ORIENTATION,FACE)) //then correct the normal
Kirill Terekhov's avatar
Kirill Terekhov committed
365
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
366
367
368
			Element::adj_type & lc = LowConn(e);
			for(Element::adj_type::iterator it = lc.begin(); it != lc.end(); ++it)
				if( !GetMarker(*it,HideMarker()) && HighConn(*it).size() == 1 )
Kirill Terekhov's avatar
Kirill Terekhov committed
369
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
370
					Face(this,*it)->FixNormalOrientation();
Kirill Terekhov's avatar
Kirill Terekhov committed
371
372
373
374
				}
		}
		for(d = MEASURE; d <= BARYCENTER; d++) // compute the rest
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
375
			if( HaveGeometricData(d,GetHandleElementType(e)) )
Kirill Terekhov's avatar
Kirill Terekhov committed
376
377
			{
				Tag t = GetGeometricTag(d);
Kirill Terekhov's avatar
Kirill Terekhov committed
378
379
				Storage::real * a = static_cast<Storage::real *>(MGetDenseLink(e,t));
				HideGeometricData(d,GetHandleElementType(e));
Kirill Terekhov's avatar
Kirill Terekhov committed
380
				GetGeometricData(e,d,a);
Kirill Terekhov's avatar
Kirill Terekhov committed
381
				ShowGeometricData(d,GetHandleElementType(e));
Kirill Terekhov's avatar
Kirill Terekhov committed
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
			}
		}
	}
	
	
	void Mesh::RemoveGeometricData(GeomParam table)
	{
		for(GeomParam::iterator it = table.begin(); it != table.end(); ++it)
		{
			if( it->first == MEASURE    ) 
			{
				if(measure_tag.isValid())    
					measure_tag    = DeleteTag(measure_tag   ,it->second);
				for(ElementType etype = EDGE; etype <= CELL; etype = etype << 1) if( etype & it->second) HideGeometricData(MEASURE,etype);
			}
			if( it->first == CENTROID   ) 
			{
				if(centroid_tag.isValid())   
					centroid_tag   = DeleteTag(centroid_tag  ,it->second);
				for(ElementType etype = EDGE; etype <= CELL; etype = etype << 1) if( etype & it->second) HideGeometricData(CENTROID,etype);
			}
			if( it->first == BARYCENTER ) 
			{
				if(barycenter_tag.isValid()) 
					barycenter_tag = DeleteTag(barycenter_tag,it->second);
				for(ElementType etype = EDGE; etype <= CELL; etype = etype << 1) if( etype & it->second) HideGeometricData(BARYCENTER,etype);
			}
			if( it->first == NORMAL     ) 
			{
				if(normal_tag.isValid())
					normal_tag     = DeleteTag(normal_tag    ,it->second);
				for(ElementType etype = FACE; etype <= CELL; etype = etype << 1) if( etype & it->second) HideGeometricData(NORMAL,etype);
			}
			if( it->first == ORIENTATION) 
				if( FACE & it->second) HideGeometricData(ORIENTATION,FACE);
		}
	}
	
	void Mesh::RestoreGeometricTags()
	{
		for(GeometricData gtype = MEASURE; gtype <= NORMAL; gtype++)
		{
			bool restore = false;
			for(ElementType etype = EDGE; etype <= CELL && !restore; etype = etype << 1)
				if( HaveGeometricData(gtype,etype) )
					restore = true;
			if( restore )
			{
				switch(gtype)
				{
				case MEASURE:       measure_tag = GetTag("GEOM_UTIL_MEASURE");    break;
				case CENTROID:     centroid_tag = GetTag("GEOM_UTIL_CENTROID");   break;
				case BARYCENTER: barycenter_tag = GetTag("GEOM_UTIL_BARYCENTER"); break;
				case NORMAL:         normal_tag = GetTag("GEOM_UTIL_NORMAL");     break;
				}
			}
		}
	}
	
	void Mesh::PrepareGeometricData(GeomParam table)
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
443
		std::sort(&*table.begin(),&*table.end());
Kirill Terekhov's avatar
Kirill Terekhov committed
444
445
446
447
448
449
		for(GeomParam::iterator it = table.begin(); it != table.end(); ++it)
		{
			GeometricData types = it->first;
			ElementType mask = it->second;
			if( types == ORIENTATION )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
450
				//std::cout << "ORIENTATION" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
451
				if( mask & FACE )
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
452
453
454
455
456
					for(integer e = 0; e < FaceLastLocalID(); ++e) 
					{
						if( isValidElement(FACE,e) )
							Face(this,ComposeHandle(FACE,e))->FixNormalOrientation();
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
457
458
459
460
				ShowGeometricData(ORIENTATION,FACE);
			}
			if( types == MEASURE )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
461
				//std::cout << "MEASURE" << std::endl;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
462
				for(ElementType etype = EDGE; etype <= CELL; etype = NextElementType(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
463
464
465
466
				{
					if( (mask & etype) && !HaveGeometricData(MEASURE,etype))
					{
						measure_tag = CreateTag("GEOM_UTIL_MEASURE",DATA_REAL,etype,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
467
468
469
470
471
						for(integer e = 0; e < LastLocalID(etype); ++e) if( isValidElement(etype,e) )
						{
							HandleType h = ComposeHandle(etype,e);
							GetGeometricData(h,MEASURE,static_cast<Storage::real *>(MGetDenseLink(h,measure_tag)));
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
472
473
474
475
476
477
						ShowGeometricData(MEASURE,etype);
					}
				}
			}
			if( types == CENTROID )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
478
				//std::cout << "CENTROID" << std::endl;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
479
				for(ElementType etype = EDGE; etype <= CELL; etype = NextElementType(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
480
481
482
483
				{
					if( (mask & etype) && !HaveGeometricData(CENTROID,etype))
					{
						centroid_tag = CreateTag("GEOM_UTIL_CENTROID",DATA_REAL,etype,NONE,GetDimensions());
Kirill Terekhov's avatar
Kirill Terekhov committed
484
						for(integer k = 0; k < LastLocalID(etype); ++k) if( isValidElement(etype,k) )
Kirill Terekhov's avatar
Kirill Terekhov committed
485
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
486
487
							HandleType h = ComposeHandle(etype,k);
							GetGeometricData(h,CENTROID,static_cast<Storage::real *>(MGetDenseLink(h,centroid_tag)));
Kirill Terekhov's avatar
Kirill Terekhov committed
488
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
489
490
491
492
493
494
						ShowGeometricData(CENTROID,etype);
					}
				}
			}
			if( types == BARYCENTER )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
495
				//std::cout << "BARYCENTER" << std::endl;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
496
				for(ElementType etype = EDGE; etype <= CELL; etype = NextElementType(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
497
498
499
500
				{
					if( (mask & etype) && !HaveGeometricData(BARYCENTER,etype))
					{
						barycenter_tag = CreateTag("GEOM_UTIL_BARYCENTER",DATA_REAL,etype,NONE,GetDimensions());
Kirill Terekhov's avatar
Kirill Terekhov committed
501
502
503
504
505
						for(integer e = 0; e < LastLocalID(etype); ++e) if( isValidElement(etype,e) )
						{
							HandleType h = ComposeHandle(etype,e);
							GetGeometricData(h,BARYCENTER,static_cast<Storage::real *>(MGetDenseLink(h,barycenter_tag)));
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
506
507
508
509
510
511
						ShowGeometricData(BARYCENTER,etype);
					}
				}	
			}
			if( types == NORMAL )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
512
				//std::cout << "NORMAL" << std::endl;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
513
				for(ElementType etype = FACE; etype <= CELL; etype = NextElementType(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
514
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
515
516
517
					if( (mask & etype) && !HaveGeometricData(NORMAL,etype))
					{
						normal_tag = CreateTag("GEOM_UTIL_NORMAL",DATA_REAL,etype,NONE,GetDimensions());
Kirill Terekhov's avatar
Kirill Terekhov committed
518
519
520
521
522
						for(integer e = 0; e < LastLocalID(etype); ++e) if( isValidElement(etype,e) )
						{
							HandleType h = ComposeHandle(etype,e);
							GetGeometricData(h,NORMAL,static_cast<Storage::real *>(MGetDenseLink(h,normal_tag)));
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
523
524
						ShowGeometricData(NORMAL,etype);
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
525
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
526
527
528
529
			}
		}
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
530
	void Mesh::GetGeometricData(HandleType e, GeometricData type, Storage::real * ret)
Kirill Terekhov's avatar
Kirill Terekhov committed
531
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
532
533
534
		assert(e != InvalidHandle());
		assert(ret != NULL);
		assert(type == MEASURE || type == CENTROID || type == BARYCENTER || type == NORMAL);
Kirill Terekhov's avatar
Kirill Terekhov committed
535
536
537
		ElementType etype = GetHandleElementType(e);
		integer edim = Element::GetGeometricDimension(GetGeometricType(e));
		integer mdim = GetDimensions();
Kirill Terekhov's avatar
Kirill Terekhov committed
538
539
540
541
542
		switch(type)
		{
			case MEASURE:
			if( HaveGeometricData(MEASURE,etype) )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
543
				*ret = static_cast<Storage::real *>(MGetDenseLink(e,measure_tag))[0];
Kirill Terekhov's avatar
Kirill Terekhov committed
544
545
546
547
548
549
550
551
552
				//~ if( isnan(*ret) || fabs(*ret) < 1e-15  ) throw -1;
			}
			else
			{
				switch(edim)
				{
					case 0: *ret = 0; break;
					case 1: //length of edge
					{
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
553
						ElementArray<Node> nodes = Element(this,e)->getNodes();
Kirill Terekhov's avatar
Kirill Terekhov committed
554
555
556
557
558
559
560
						if( nodes.size() > 1 )
						{
							Storage::real c[3];
							vec_diff(nodes[0]->Coords().data(),nodes[1]->Coords().data(),c,mdim);
							*ret = vec_len(c,mdim);
						}
						else *ret = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
561
562
563
564
565
						//~ if( isnan(*ret) || fabs(*ret) < 1e-15  ) throw -1;
						break;
					}
					case 2: //area of face
					{
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
566
						ElementArray<Node> nodes = Element(this,e)->getNodes();
Kirill Terekhov's avatar
Kirill Terekhov committed
567
						if( nodes.size() > 2 )
Kirill Terekhov's avatar
Kirill Terekhov committed
568
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
569
570
571
							real x[3] = {0,0,0};
							Storage::real_array x0 = nodes[0].Coords();
							for(ElementArray<Node>::size_type i = 1; i < nodes.size()-1; i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
572
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
573
574
575
576
577
578
579
580
								Storage::real_array v1 = nodes[i].Coords();
								Storage::real_array v2 = nodes[i+1].Coords();
								if( mdim == 3 )
								{
									x[0] += (v1[1]-x0[1])*(v2[2]-x0[2]) - (v1[2]-x0[2])*(v2[1]-x0[1]);
									x[1] += (v1[2]-x0[2])*(v2[0]-x0[0]) - (v1[0]-x0[0])*(v2[2]-x0[2]);
								}
								x[2] += (v1[0]-x0[0])*(v2[1]-x0[1]) - (v1[1]-x0[1])*(v2[0]-x0[0]);
Kirill Terekhov's avatar
Kirill Terekhov committed
581
							}
Kirill Terekhov's avatar
Kirill Terekhov committed
582
583
							*ret = sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])*0.5;
						} else *ret = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
584
585
586
587
588
						//~ if( isnan(*ret) || fabs(*ret) < 1e-15  ) throw -1;
						break;
					}
					case 3: //volume of cell
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
589
590
						Cell me = Cell(this,e);
						ElementArray<Face> faces = me->getFaces();
Kirill Terekhov's avatar
Kirill Terekhov committed
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
						*ret = 0;
						/*
						Storage::real d;
						for(unsigned i = 0; i < faces.size(); i++)
						{
							d = 0;
							adjacent<Node> nodes = faces[i].getNodes();
							Storage::real_array a = nodes[0].Coords();
							for(unsigned j = 1; j < nodes.size()-1; j++)
							{
								Storage::real_array b = nodes[j].Coords();
								Storage::real_array c = nodes[j+1].Coords();
								d += det3v(&a[0],&b[0],&c[0]);
							}
							*ret += (2*faces[i].FaceOrientedOutside(e->getAsCell())-1)*d;
						}
						*ret /= 6.0;
						*/
Kirill Terekhov's avatar
Kirill Terekhov committed
609
						if( faces.size() > 3 )
Kirill Terekhov's avatar
Kirill Terekhov committed
610
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
611
612
							Storage::real fcnt[3], fnrm[3];// , area;
							for(ElementArray<Face>::size_type i = 0; i < faces.size(); i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
613
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
614
615
616
617
618
619
620
621
622
623
624
625
626
								faces[i]->Centroid(fcnt);
								faces[i]->OrientedNormal(me,fnrm);
								/*
								area = sqrt(vec_dot_product(fnrm,fnrm,mdim));
								if( area > 0 )
								{
									fnrm[0] /= area;
									fnrm[1] /= area;
									fnrm[2] /= area;
									*ret += vec_dot_product(fcnt,fnrm,mdim) * area;
								}
								*/
								*ret += vec_dot_product(fcnt,fnrm,mdim);
Kirill Terekhov's avatar
Kirill Terekhov committed
627
							}
Kirill Terekhov's avatar
Kirill Terekhov committed
628
							*ret /= 3.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
						}
						/*
						if( *ret < 0  || (*ret) != (*ret) )
						{
							Storage::real test = 0;
							
							
							//~ Storage::real fcnt[3], fnrm[3], area;
							for(unsigned i = 0; i < faces.size(); i++)
							{
								if( faces[i].FixNormalOrientation() ) std::cout << faces[i].LocalID() << " normal refixed " << faces[i].nbAdjElements(CELL) << std::endl;
								//~ faces[i].Centroid(fcnt);
								//~ faces[i].OrientedNormal(e->getAsCell(),fnrm);
								//~ area = sqrt(vec_dot_product(fnrm,fnrm,mdim));
								//~ fnrm[0] /= area;
								//~ fnrm[1] /= area;
								//~ fnrm[2] /= area;
								//~ test += vec_dot_product(fcnt,fnrm,mdim) * area / 3.0;
							}
							
							//~ e->Centroid(fcnt);
							
							Storage::real d;
							for(unsigned i = 0; i < faces.size(); i++)
							{
								d = 0;
								adjacent<Node> nodes = faces[i].getNodes();
								Storage::real_array a = nodes[0].Coords();
								for(unsigned j = 1; j < nodes.size()-1; j++)
								{
									Storage::real_array b = nodes[j].Coords();
									Storage::real_array c = nodes[j+1].Coords();
									d += det3v(&a[0],&b[0],&c[0]);
								}
								test += (2*faces[i].FaceOrientedOutside(e->getAsCell())-1)*d;
							}
							test /= 6.0;
						
						
							std::cout << "alg1 " << *ret << " alg2 " << test << " on " << Element::GeometricTypeName(e->GetGeometricType()) << " " << e->GlobalID() << " " << e->LocalID() << std::endl;
							//e->Integer(tag_topologyerror) = 1;
						}
						*/
						//~ if( *ret < 0 )
						//~ {
							//~ std::cout << "negative volume! " << *ret << std::endl;
							//~ std::cout << "element " << Element::GeometricTypeName(e->GetGeometricType()) << " faces " << faces.size() << std::endl;
							//~ *ret = 0;
							//~ for(unsigned i = 0; i < faces.size(); i++)
							//~ {
								//~ std::cout << "face " << i << "/" << faces.size() << std::endl;
								//~ d = 0;
								//~ adjacent<Node> nodes = faces[i].getNodes();
								//~ Storage::real_array a = nodes[0].Coords();
								//~ std::cout << "node 0 " << a[0] << " " << a[1] << " " << a[2] << " id " << nodes[0].LocalID() << std::endl; 
								//~ for(unsigned j = 1; j < nodes.size()-1; j++)
								//~ {
									//~ Storage::real_array b = nodes[j].Coords();
									//~ std::cout << "node " << j << " " << b[0] << " " << b[1] << " " << b[2] << " id " << nodes[j].LocalID() << std::endl;
									//~ Storage::real_array c = nodes[j+1].Coords();
									//~ d += det3v(&a[0],&b[0],&c[0]);
								//~ }
								//~ a = nodes[nodes.size()-1].Coords();
								//~ std::cout << "node " << nodes.size()-1 << " " << a[0] << " " << a[1] << " " << a[2] << " id " << nodes[nodes.size()-1].LocalID() << std::endl; 
								//~ std::cout << "d = " << d << std::endl;
								//~ std::cout << "orientation = " << (2*faces[i].FaceOrientedOutside(e->getAsCell())-1) << std::endl;
								//~ Storage::real old = *ret, add = (2*faces[i].FaceOrientedOutside(e->getAsCell())-1)*d;
								//~ *ret = old + add;
								//~ std::cout << old << " + " << add << " = " << *ret << std::endl;
							//~ }
							//~ std::cout << "result " << *ret/6.0 << std::endl;
							//~ std::cout << "trying with fix " << std::endl;
							//~ *ret = 0;
							//~ for(unsigned i = 0; i < faces.size(); i++)
							//~ {
								//~ std::cout << "face " << i << "/" << faces.size() << " fixed " << faces[i].FixNormalOrientation() << " cells " << faces[i].nbAdjElements(CELL) << std::endl;
								//~ d = 0;
								//~ adjacent<Node> nodes = faces[i].getNodes();
								//~ Storage::real_array a = nodes[0].Coords();
								//~ std::cout << "node 0 " << a[0] << " " << a[1] << " " << a[2] << " id " << nodes[0].LocalID() << std::endl; 
								//~ for(unsigned j = 1; j < nodes.size()-1; j++)
								//~ {
									//~ Storage::real_array b = nodes[j].Coords();
									//~ std::cout << "node " << j << " " << b[0] << " " << b[1] << " " << b[2] << " id " << nodes[j].LocalID() << std::endl;
									//~ Storage::real_array c = nodes[j+1].Coords();
									//~ d += det3v(&a[0],&b[0],&c[0]);
								//~ }
								//~ a = nodes[nodes.size()-1].Coords();
								//~ std::cout << "node " << nodes.size()-1 << " " << a[0] << " " << a[1] << " " << a[2] << " id " << nodes[nodes.size()-1].LocalID() << std::endl; 
								//~ std::cout << "d = " << d << std::endl;
								//~ std::cout << "orientation = " << (2*faces[i].FaceOrientedOutside(e->getAsCell())-1) << std::endl;
								//~ Storage::real old = *ret, add = (2*faces[i].FaceOrientedOutside(e->getAsCell())-1)*d;
								//~ *ret = old + add;
								//~ std::cout << old << " + " << add << " = " << *ret << std::endl;
							//~ }
							//~ std::cout << "result " << *ret/6.0 << std::endl;
						//~ }
						
						//~ if( isnan(*ret) || fabs(*ret) < 1e-15  ) throw -1;
						break;
					}
				}
				//~ throw -1;
				//~ if( (*ret) != (*ret) || *ret < 0  ) 
				//~ {
					//~ std::cout << "bad measure: " << *ret << " for " << ElementTypeName(e->GetElementType()) << " " << Element::GeometricTypeName(e->GetGeometricType()) << " edim " << edim << std::endl;
					//~ 
				//~ }
				
			}
			//~ if( isnan(*ret) || fabs(*ret) < 1e-15  ) throw -1;
			break;
			case CENTROID:
			if(etype == NODE )
Kirill Terekhov's avatar
Kirill Terekhov committed
743
				memcpy(ret,MGetDenseLink(e,CoordsTag()),sizeof(real)*mdim);
Kirill Terekhov's avatar
Kirill Terekhov committed
744
745
			else if(HaveGeometricData(CENTROID,etype))
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
746
				memcpy(ret,MGetDenseLink(e,centroid_tag),sizeof(real)*mdim);
Kirill Terekhov's avatar
Kirill Terekhov committed
747
748
749
			}
			else
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
750
				ElementArray<Node> nodes = Element(this,e)->getNodes();
Kirill Terekhov's avatar
Kirill Terekhov committed
751
				memset(ret,0,sizeof(real)*mdim);
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
752
				if(nodes.size() != 0)
Kirill Terekhov's avatar
Kirill Terekhov committed
753
				{
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
754
755
756
757
758
759
760
					Storage::real div = 1.0/nodes.size();
					for(ElementArray<Node>::size_type i = 0; i < nodes.size(); i++)
					{
						Storage::real_array c =nodes[i].Coords();
						for(integer j = 0; j < mdim; j++) ret[j] += c[j];
					}
					for(integer j = 0; j < mdim; j++) ret[j] *= div;
Kirill Terekhov's avatar
Kirill Terekhov committed
761
762
763
764
765
				}
			}
			break;
			case BARYCENTER:
			if( etype == NODE )
Kirill Terekhov's avatar
Kirill Terekhov committed
766
				memcpy(ret,MGetDenseLink(e,CoordsTag()),sizeof(real)*mdim);
Kirill Terekhov's avatar
Kirill Terekhov committed
767
			else if(HaveGeometricData(BARYCENTER,etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
768
				memcpy(ret,MGetDenseLink(e,barycenter_tag),sizeof(real)*mdim);
Kirill Terekhov's avatar
Kirill Terekhov committed
769
770
771
772
773
			else
			{
				memset(ret,0,sizeof(real)*mdim);
				if( edim == 1 )
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
774
					ElementArray<Node> n = Element(this,e)->getNodes();
Kirill Terekhov's avatar
Kirill Terekhov committed
775
776
777
778
779
780
781
782
783
784
785
786
					if( n.size() == 2 )
					{
						Storage::real_array a = n[0].Coords();
						Storage::real_array b = n[1].Coords();
						for(integer j = 0; j < dim; j++) 
							ret[j] = (a[j] + b[j])*0.5;
					}
					else if( n.size() == 1 )
					{
						Storage::real_array a = n[0].Coords();
						for(integer j = 0; j < dim; j++) ret[j] = a[j];
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
787
788
789
				}
				else if( edim == 2 )
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
790
					ElementArray<Node> nodes = Element(this,e)->getNodes();
Kirill Terekhov's avatar
Kirill Terekhov committed
791
792
793
794
795
796
797
798
799
					real s,d, x1[3] = {0,0,0},x2[3] = {0,0,0},x[3] = {0,0,0};
					//here we compute area of polygon
					//~ if( HaveGeometricData(MEASURE,etype) && HaveGeometricData(NORMAL,etype) )
					//~ {
						//~ s = e->RealDF(measure_tag);
						//~ memcpy(x,&e->RealDF(normal_tag),sizeof(Storage::real)*mdim);
					//~ }
					//~ else
					//~ {
Kirill Terekhov's avatar
Kirill Terekhov committed
800
					if( nodes.size() > 2 )
Kirill Terekhov's avatar
Kirill Terekhov committed
801
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
802
803
						Storage::real_array x0 = nodes[0].Coords();
						for(ElementArray<Node>::size_type i = 1; i < nodes.size()-1; i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
804
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
805
806
807
808
809
810
811
812
							Storage::real_array v1 = nodes[i].Coords();
							Storage::real_array v2 = nodes[i+1].Coords();
							if( mdim == 3 )
							{
								x[0] += (v1[1]-x0[1])*(v2[2]-x0[2]) - (v1[2]-x0[2])*(v2[1]-x0[1]);
								x[1] += (v1[2]-x0[2])*(v2[0]-x0[0]) - (v1[0]-x0[0])*(v2[2]-x0[2]);
							}
							x[2] += (v1[0]-x0[0])*(v2[1]-x0[1]) - (v1[1]-x0[1])*(v2[0]-x0[0]);
Kirill Terekhov's avatar
Kirill Terekhov committed
813
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
						s = sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
						x[0] /= s; x[1] /= s; x[2] /= s; //here we obtain the unit normal
						//~ }
						//here we compute the center
						Storage::real_array v0 = nodes[0].Coords();
						for(ElementArray<Node>::size_type j = 1; j < nodes.size()-1; j++)
						{
							Storage::real_array v1 = nodes[j].Coords();
							Storage::real_array v2 = nodes[j+1].Coords();
							for(integer k = 0; k < mdim; k++)
							{
								x1[k] = v0[k] - v1[k];
								x2[k] = v0[k] - v2[k];
							}
							d = det3v(x1,x2,x); //here we use unit normal
							for(integer k = 0; k < mdim; k++) 
								ret[k] += d*(v0[k]+v1[k]+v2[k]);
						}
						for(integer k = 0; k < mdim; k++) ret[k] /= 3.0 * s;
Kirill Terekhov's avatar
Kirill Terekhov committed
833
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
834
					else
Kirill Terekhov's avatar
Kirill Terekhov committed
835
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
836
						for(ElementArray<Node>::size_type i = 0; i < nodes.size(); i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
837
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
838
839
							Storage::real_array c = nodes[i].Coords();
							for(integer k = 0; k < mdim; k++) ret[k] += c[k];
Kirill Terekhov's avatar
Kirill Terekhov committed
840
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
841
						for(integer k = 0; k < mdim; k++) ret[k] /= static_cast<Storage::real>(nodes.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
842
843
844
845
					}
				}
				else if( edim == 3 )
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
846
					ElementArray<Face> faces = Element(this,e)->getFaces();
Kirill Terekhov's avatar
Kirill Terekhov committed
847
					real d,c,vol = 0, y[3];
Kirill Terekhov's avatar
Kirill Terekhov committed
848
					for(ElementArray<Face>::size_type i = 0; i < faces.size(); i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
849
850
					{
						d = y[0] = y[1] = y[2] = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
851
						ElementArray<Node> nodes = faces[i].getNodes();
Kirill Terekhov's avatar
Kirill Terekhov committed
852
						Storage::real_array v0 = nodes[0].Coords();
Kirill Terekhov's avatar
Kirill Terekhov committed
853
						for(ElementArray<Node>::size_type j = 1; j < nodes.size()-1; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
854
855
856
						{
							Storage::real_array v1 = nodes[j].Coords();
							Storage::real_array v2 = nodes[j+1].Coords();
Kirill Terekhov's avatar
Kirill Terekhov committed
857
							c = det3v(v0.data(),v1.data(),v2.data());
Kirill Terekhov's avatar
Kirill Terekhov committed
858
859
860
861
862
							d += c;
							y[0] += c * (v0[0] + v1[0] + v2[0]);
							y[1] += c * (v0[1] + v1[1] + v2[1]);
							y[2] += c * (v0[2] + v1[2] + v2[2]);
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
863
						c = faces[i].FaceOrientedOutside(Cell(this,e)) ? 1 : -1;
Kirill Terekhov's avatar
Kirill Terekhov committed
864
865
866
867
868
869
870
871
872
873
874
875
876
						ret[0] += c * y[0];
						ret[1] += c * y[1];
						ret[2] += c * y[2];
						vol += c*d;
					}
					ret[0] /= vol*4;
					ret[1] /= vol*4;
					ret[2] /= vol*4;
				}
			}
			break;
			case NORMAL:
			if( HaveGeometricData(NORMAL,etype) )
Kirill Terekhov's avatar
Kirill Terekhov committed
877
				memcpy(ret,MGetDenseLink(e,normal_tag),sizeof(real)*mdim);
Kirill Terekhov's avatar
Kirill Terekhov committed
878
879
880
			else
			{
				memset(ret,0,sizeof(real)*mdim);
Kirill Terekhov's avatar
Kirill Terekhov committed
881
				if( edim == 2 )//&& mdim == 3)
Kirill Terekhov's avatar
Kirill Terekhov committed
882
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
883
					ElementArray<Node> nodes = Element(this,e)->getNodes();
Kirill Terekhov's avatar
Kirill Terekhov committed
884
885
					
					Storage::real_array x0 = nodes[0].Coords();
Kirill Terekhov's avatar
Kirill Terekhov committed
886
					for(ElementArray<Node>::size_type i = 1; i < nodes.size()-1; i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
					{
						Storage::real_array a = nodes[i].Coords();
						Storage::real_array b = nodes[i+1].Coords();
						ret[0] += (a[1]-x0[1])*(b[2]-x0[2]) - (a[2]-x0[2])*(b[1]-x0[1]);
						ret[1] += (a[2]-x0[2])*(b[0]-x0[0]) - (a[0]-x0[0])*(b[2]-x0[2]);
						ret[2] += (a[0]-x0[0])*(b[1]-x0[1]) - (a[1]-x0[1])*(b[0]-x0[0]);
					}
					/*
					for(unsigned i = 0; i < nodes.size(); i++)
					{
						Storage::real_array a = nodes[i].Coords();
						Storage::real_array b = nodes[(i+1)%nodes.size()].Coords();
						ret[0] += a[1]*b[2] - a[2]*b[1];
						ret[1] += a[2]*b[0] - a[0]*b[2];
						ret[2] += a[0]*b[1] - a[1]*b[0];
					}
					*/
					ret[0] *= 0.5;
					ret[1] *= 0.5;
					ret[2] *= 0.5;
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
908
				else if( edim == 1 )//&& mdim == 2 )
Kirill Terekhov's avatar
Kirill Terekhov committed
909
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
910
					ElementArray<Node> nodes = Element(this,e)->getNodes();
Kirill Terekhov's avatar
Kirill Terekhov committed
911
					if( nodes.size() > 1 )
Kirill Terekhov's avatar
Kirill Terekhov committed
912
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
913
914
915
916
917
918
919
920
921
922
923
924
925
						Storage::real_array a = nodes[0].Coords();
						Storage::real_array b = nodes[1].Coords();
						ret[0] = b[1] - a[1];
						ret[1] = a[0] - b[0];
						Storage::real l = sqrt(ret[0]*ret[0]+ret[1]*ret[1]);
						if( l )
						{
							ret[0] /= l;
							ret[1] /= l;
						}
						l = sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1]));
						ret[0] *= l;
						ret[1] *= l;
Kirill Terekhov's avatar
Kirill Terekhov committed
926
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
927
928
929
930
931
932
933
934
935
936
937
938
				}
			}
			break;
		}
		//~ if( type == MEASURE )
		//~ {
			//~ if( isnan(*ret) || fabs(*ret) < 1e-15  ) throw -1;
		//~ }
	}
	


Kirill Terekhov's avatar
Kirill Terekhov committed
939
	bool Element::Planarity() const
Kirill Terekhov's avatar
Kirill Terekhov committed
940
941
	{
		Mesh * m = GetMeshLink();
Kirill Terekhov's avatar
Kirill Terekhov committed
942
		integer dim = m->GetDimensions();
Kirill Terekhov's avatar
Kirill Terekhov committed
943
		if( dim < 3 ) return true;
Kirill Terekhov's avatar
Kirill Terekhov committed
944
		ElementArray<Node> p = getNodes();
Kirill Terekhov's avatar
Kirill Terekhov committed
945
		if( p.size() <= 3 ) return true;
Kirill Terekhov's avatar
Kirill Terekhov committed
946
		ElementArray<Node>::size_type i, s = p.size();
Kirill Terekhov's avatar
Kirill Terekhov committed
947
		Storage::real v[2][3] = {{0,0,0},{0,0,0}};
Kirill Terekhov's avatar
Kirill Terekhov committed
948
949
		vec_diff(p[1].Coords().data(),p[0].Coords().data(),v[0],3);
		vec_diff(p[2].Coords().data(),p[0].Coords().data(),v[1],3);
Kirill Terekhov's avatar
Kirill Terekhov committed
950
951
952
		vec_cross_product(v[0],v[1],v[1]);
		for(i = 3; i < s; i++)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
953
			vec_diff(p[i].Coords().data(),p[0].Coords().data(),v[0],3);
Kirill Terekhov's avatar
Kirill Terekhov committed
954
955
956
957
958
959
			if( fabs(vec_dot_product(v[0],v[1],3)) > m->GetEpsilon() ) return false;
		}
		return true;
	}


Kirill Terekhov's avatar
Kirill Terekhov committed
960
	bool Cell::Closure() const
Kirill Terekhov's avatar
Kirill Terekhov committed
961
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
962
963
		adj_type & lc = GetMeshLink()->LowConn(GetHandle());
		return lc.size() > 0 ? GetMeshLink()->TestClosure(lc.data(),static_cast<integer>(lc.size())) : false;
Kirill Terekhov's avatar
Kirill Terekhov committed
964
965
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
966
	bool Face::Closure() const
Kirill Terekhov's avatar
Kirill Terekhov committed
967
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
968
969
		adj_type & lc = GetMeshLink()->LowConn(GetHandle());
		return lc.size() > 0 ? GetMeshLink()->TestClosure(lc.data(),static_cast<integer>(lc.size())) : false;
Kirill Terekhov's avatar
Kirill Terekhov committed
970
971
972
973
974
975
976
977
	}



	

	

Kirill Terekhov's avatar
Kirill Terekhov committed
978
	bool Element::Boundary() const
Kirill Terekhov's avatar
Kirill Terekhov committed
979
980
981
982
983
984
985
986
987
988
	{
		switch(GetElementType())
		{
			case FACE:
				if( nbAdjElements(CELL) == 1 )
				{
					if( getAsFace()->BackCell()->GetStatus() != Element::Ghost )
						return true;
				}
				return false;
Kirill Terekhov's avatar
Kirill Terekhov committed
989
			case CELL:
Kirill Terekhov's avatar
Kirill Terekhov committed
990
991
992
			case EDGE:
			case NODE:
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
993
994
				ElementArray<Element> faces = getAdjElements(FACE);
				for(ElementArray<Element>::iterator it = faces.begin(); it != faces.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
995
996
997
998
999
1000
1001
1002
1003
					if( it->Boundary() ) return true;
				return false;
			}
			default: return false;
		}
		return false;
	}
	
	
Kirill Terekhov's avatar
Kirill Terekhov committed
1004
	bool Face::CheckNormalOrientation() const
Kirill Terekhov's avatar
Kirill Terekhov committed
1005
1006
	{
		Mesh * m = GetMeshLink();
Kirill Terekhov's avatar
Kirill Terekhov committed
1007
1008
1009
1010
		integer dim = m->GetDimensions();
		Cell c1 = BackCell();
		Cell c2 = FrontCell();
		if( c1.isValid() )
Kirill Terekhov's avatar
Kirill Terekhov committed
1011
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1012
			real nf[3], cf[3],fc[3], cc[3], dot, dot2;
Kirill Terekhov's avatar
Kirill Terekhov committed
1013
1014
1015
1016
1017
1018
1019
1020
			UnitNormal(nf);
			Centroid(fc);
			c1->Centroid(cc);
			//~ if( !c->Inside(cc) ) std::cout << "centroid is outside of " << Element::GeometricTypeName(c->GetGeometricType()) << std::endl;
			vec_diff(cc,fc,cf,dim);
			vec_normalize(cf,dim);
			dot = vec_dot_product(nf,cf,dim);
			//~ c = FrontCell();
Kirill Terekhov's avatar
Kirill Terekhov committed
1021
			if( fabs(dot) < 0.25 && c2.isValid() )
Kirill Terekhov's avatar
Kirill Terekhov committed
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
			{
				c2->Centroid(cc);
				vec_diff(cc,fc,cf,dim);
				vec_normalize(cf,dim);
				dot2 = vec_dot_product(nf,cf,dim);
				if( fabs(dot2) > fabs(dot) )
					dot = -dot2;
			}
			if( dot > 0 )
				return false;
		}
		return true;
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
1036
	bool Face::FixNormalOrientation() const
Kirill Terekhov's avatar
Kirill Terekhov committed
1037
1038
1039
1040
1041
1042
	{
		if( !CheckNormalOrientation() )
		{	
			ReorderEdges();
			if( GetMeshLink()->HaveGeometricData(NORMAL,FACE) )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1043
1044
1045
				real_array nrm = GetMeshLink()->RealArrayDF(GetHandle(),GetMeshLink()->GetGeometricTag(NORMAL));
				for(real_array::size_type it = 0; it < nrm.size(); ++it) 
					nrm[it] = -nrm[it];
Kirill Terekhov's avatar
Kirill Terekhov committed
1046
1047
1048
1049
1050
1051
			}
			return true;
		}
		return false;
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
1052
	Storage::real meantri(Storage::real * v0, Storage::real * v1, Storage::real * v2, Storage::integer dim, Storage::real (*func)(Storage::real* x,Storage::real), Storage::real time)
Kirill Terekhov's avatar
Kirill Terekhov committed
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
	{
		Storage::real value = 0;
		static const Storage::real w[4] =   { -0.149570044467670, 0.175615257433204, 0.053347235608839 , 0.077113760890257};
		static const Storage::real a[4][3] =
		{
			{0.333333333333333,0.333333333333333,0.333333333333333},
			{0.479308067841923,0.260345966079038,0.260345966079038},
			{0.869739794195568,0.065130102902216,0.065130102902216},
			{0.638444188569809,0.312865496004875,0.048690315425316}
		};
		Storage::real XYG[13][3];
Kirill Terekhov's avatar
Kirill Terekhov committed
1064
		for (Storage::integer i = 0 ; i < dim; i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1065
1066
1067
1068
			XYG[0][i] = 0.33333333333333333333*(v0[i]+v1[i]+v2[i]);
		 value += w[0] * func(XYG[0],time);
		for (int i = 0 ; i < 3 ; i++ )
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1069
			for (Storage::integer j = 0 ; j < dim; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1070
1071
1072
1073
1074
				XYG[1+i][j] = v0[j] + (v1[j] - v0[j]) * a[1][i] + (v2[j] - v0[j])*a[1][(i+1)%3];
			value += w[1] * func(XYG[1+i],time);
		}
		for (int i = 0 ; i < 3 ; i++ )
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1075
			for (Storage::integer j = 0 ; j < dim; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1076
1077
1078
1079
1080
				XYG[4+i][j] = v0[j] + (v1[j] - v0[j]) * a[2][i] + (v2[j] - v0[j])*a[2][(i+1)%3];
			value += w[2] * func(XYG[4+i],time);
		}
		for (int i = 0 ; i < 3 ; i++ )
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1081
			for (Storage::integer j = 0 ; j < dim; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
			{
				XYG[7+2*i][j] = v0[j] + (v1[j] - v0[j]) * a[3][i] + (v2[j] - v0[j])*a[3][(i+1)%3];
				XYG[8+2*i][j] = v0[j] + (v1[j] - v0[j]) * a[3][(i+1)%3] + (v2[j] - v0[j])*a[3][i];
			}
			value += w[3] * func(XYG[7+2*i],time);
			value += w[3] * func(XYG[8+2*i],time);
		}
		return value;
	}

	Storage::real meantet(Storage::real * v0, Storage::real * v1, Storage::real * v2, Storage::real * v3, Storage::real (*func)(Storage::real* x,Storage::real),Storage::real time)
	{
		Storage::real value;
		static const Storage::real T5A = 0.25, W5A = 0.11851851851852;
		static const Storage::real T5B = 0.09197107805272, T5C = 0.72408676584183, W5B = 0.07193708377902;
		static const Storage::real T5D = 0.31979362782963, T5E = 0.04061911651111;
		static const Storage::real W5C = 0.06906820722627;
		static const Storage::real T5F = 0.05635083268963, T5G = 0.44364916731037;
		static const Storage::real W5D = 0.05291005291005;
		static const Storage::real w[15] = {W5A,W5B,W5B,W5B,W5B,W5C,W5C,W5C,W5C,W5D,W5D,W5D,W5D,W5D,W5D};
		Storage::real XYG[15][3];
		for (int i = 0; i < 3 ;i++)
		{
			XYG[0][i] = T5A * (v0[i] + v1[i] + v2[i] + v3[i]);
			XYG[1][i] = T5C * v0[i] + T5B * (v1[i]+v2[i]+v3[i]);
			XYG[2][i] = T5C * v1[i] + T5B * (v0[i]+v2[i]+v3[i]);
			XYG[3][i] = T5C * v2[i] + T5B * (v0[i]+v1[i]+v3[i]);
			XYG[4][i] = T5C * v3[i] + T5B * (v0[i]+v1[i]+v2[i]);

			XYG[5][i] = T5E * v0[i] + T5D * (v1[i]+v2[i]+v3[i]);
			XYG[6][i] = T5E * v1[i] + T5D * (v0[i]+v2[i]+v3[i]);
			XYG[7][i] = T5E * v2[i] + T5D * (v0[i]+v1[i]+v3[i]);
			XYG[8][i] = T5E * v3[i] + T5D * (v0[i]+v1[i]+v2[i]);

			XYG[9][i] = T5F * (v0[i] + v1[i]) + T5G * (v2[i] + v3[i]);
			XYG[10][i] = T5G * (v0[i] + v1[i]) + T5F * (v2[i] + v3[i]);
			XYG[11][i] = T5F * (v0[i] + v3[i]) + T5G * (v1[i] + v2[i]);
			XYG[12][i] = T5G * (v0[i] + v3[i]) + T5F * (v1[i] + v2[i]);
			XYG[13][i] = T5F * (v0[i] + v2[i]) + T5G * (v1[i] + v3[i]);
			XYG[14][i] = T5G * (v0[i] + v2[i]) + T5F * (v1[i] + v3[i]);
		}
		value = 0;
		for (int i = 0 ; i < 15 ; i++)
			value += w[i] * func(XYG[i],time);
		return value;
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
1129
	Storage::real Element::Mean(Storage::real (*func)(Storage::real* x,Storage::real),Storage::real time) const
Kirill Terekhov's avatar
Kirill Terekhov committed
1130
1131
1132
1133
	{
		Mesh * m = GetMeshLink();
		if( GetElementDimension() == 2 )
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1134
1135
1136
1137
1138
1139
1140
1141
			integer dim = m->GetDimensions();
			real val = 0, vol = 0, tvol,tval;
			real normal[3];
			real v1[3] = {0,0,0},v2[3] = {0,0,0}, product[3] = {0,0,0};
			m->GetGeometricData(GetHandle(),NORMAL,normal);
			ElementArray<Node> nodes = getNodes();
			real_array av0 = nodes.front().Coords();
			for(ElementArray<Node>::iterator it = ++nodes.begin(); it != nodes.end(); it++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1142
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1143
				ElementArray<Node>::iterator jt = it++;
Kirill Terekhov's avatar
Kirill Terekhov committed
1144
				if( it == nodes.end() ) break;
Kirill Terekhov's avatar
Kirill Terekhov committed
1145
1146
1147
				real_array av1 = jt->Coords();
				real_array av2 = it->Coords();
				tval = meantri(av0.data(),av1.data(),av2.data(),dim,func,time);
Kirill Terekhov's avatar
Kirill Terekhov committed
1148
1149
1150
				vec_diff(av1.data(),av0.data(),v1,dim);
				vec_diff(av2.data(),av0.data(),v2,dim);
				vec_cross_product(v1,v2,product);
Kirill Terekhov's avatar
Kirill Terekhov committed
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
				tvol = vec_dot_product(product,normal,dim)*0.5;
				val += tval*tvol;
				vol += tvol;
				it = jt;
			}
			return val / vol;
		}
		else if( GetElementDimension() == 3 )
		{
			
Kirill Terekhov's avatar
Kirill Terekhov committed
1161
1162
1163
			ElementArray<Element> rfaces = getAdjElements(FACE);
			array<int> n(static_cast<array<int>::size_type>(rfaces.size()));
			array<real> v;
Kirill Terekhov's avatar
Kirill Terekhov committed
1164
			int k = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1165
			for(ElementArray<Element>::iterator f = rfaces.begin(); f != rfaces.end(); f++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1166
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1167
1168
				ElementArray<Node> nodes = f->getNodes();
				int nn = n[k] = static_cast<int>(nodes.size());
Kirill Nikitin's avatar
Kirill Nikitin committed
1169
				for(int i = 0; i < nn; i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1170
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
1171
					real_array a = nodes[i].Coords();
Kirill Terekhov's avatar
Kirill Terekhov committed
1172
1173
1174
1175
1176
					v.insert(v.end(),a.begin(),a.end());
				}
				k++;
			}
			int j = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1177
1178
			real x[3] = {0,0,0}, y[3], d, vol = 0, c;
			for(array<int>::size_type i = 0; i < n.size(); i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1179
1180
			{
				y[0] = y[1] = y[2] = d = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1181
				for(array<int>::size_type k = 1; k < static_cast<array<int>::size_type>(n[i] - 1); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1182
1183
1184
1185
1186
1187
				{
					d += c = det3v(&v[j*3],&v[(j+k)*3],&v[(j+k+1)*3]);
					y[0] += c * (v[j*3+0] + v[(j+k)*3+0] + v[(j+k+1)*3+0]);
					y[1] += c * (v[j*3+1] + v[(j+k)*3+1] + v[(j+k+1)*3+1]);
					y[2] += c * (v[j*3+2] + v[(j+k)*3+2] + v[(j+k+1)*3+2]);
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
1188
				int orient = rfaces[static_cast<ElementArray<Element>::size_type>(i)]->getAsFace()->FaceOrientedOutside(getAsCell())?1:-1;
Kirill Terekhov's avatar
Kirill Terekhov committed
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
				x[0] += orient*y[0];
				x[1] += orient*y[1];
				x[2] += orient*y[2];
				vol += orient*d;
				j += n[i];
			}
			x[0] /= 4.0 * vol;
			x[1] /= 4.0 * vol;
			x[2] /= 4.0 * vol;
			vol /= 6;
			j = 0;
			vol = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1201
1202
1203
			real tvol, tval, val = 0;
			real vv0[3], vv1[3], vv2[3], prod[3];
			for(array<int>::size_type i = 0; i < n.size(); i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1204
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1205
				for(array<int>::size_type k = 1; k < static_cast<array<int>::size_type>(n[i] - 1); k++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1206
1207
1208
1209
1210
				{
					vec_diff(&v[j*3],x,vv0,3);
					vec_diff(&v[(j+k)*3],x,vv1,3);
					vec_diff(&v[(j+k+1)*3],x,vv2,3);
					vec_cross_product(vv1,vv2,prod);
Kirill Terekhov's avatar
Kirill Terekhov committed
1211
					tvol = vec_dot_product(vv0,prod,3)/6.0 * (rfaces[static_cast<ElementArray<Element>::size_type>(i)]->getAsFace()->FaceOrientedOutside(getAsCell())?1:-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
					tval = meantet(x,&v[j*3],&v[(j+k)*3],&v[(j+k+1)*3],func,time);
					val += tval * tvol;
					vol += tvol;
				}
				j+=n[i];
			}
			return val / vol;
		}
		else if ( GetElementDimension() == 1 ) //Mean value over line.
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1222
1223
1224
1225
1226
1227
			ElementArray<Node> nodes = getNodes();
			integer dim = m->GetDimensions();
			real_array x1 = nodes[0].Coords();
			real_array x2 = nodes[1].Coords();
			real middle[3];
			for (integer i = 0 ; i < dim ; i++) middle[i] = (x1[i]+x2[i])*0.5;
Kirill Terekhov's avatar
Kirill Terekhov committed
1228
			//Simpson formula
Kirill Terekhov's avatar
Kirill Terekhov committed
1229
			return (func(x1.data(),time) + 4*func(middle,time) + func(x2.data(),time))/6.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1230
1231
1232
1233
		}
		return 0;
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
	const Tag & Mesh::GetGeometricTag(GeometricData type) const
	{
		switch(type) 
		{
		case    MEASURE: return    measure_tag; 
		case   CENTROID: return   centroid_tag; 
		case BARYCENTER: return barycenter_tag; 
		case     NORMAL: return     normal_tag;
		}
		assert(false);
		static const Tag t; //invalid tag
		return t;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
1247
	
Kirill Terekhov's avatar
Kirill Terekhov committed
1248
	ElementArray<Face> Mesh::GatherBoundaryFaces()
Kirill Terekhov's avatar
Kirill Terekhov committed
1249
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
1250
1251
1252
1253
1254
1255
		ElementArray<Face> ret(this);
		for(integer f = 0; f < FaceLastLocalID(); f++) if( isValidElement(FACE,f) )
		{
			Face ef = Face(this,ComposeHandle(FACE,f));
			if( ef->Boundary() ) ret.push_back(ef);
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
1256
1257
1258
		return ret;
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
1259
	ElementArray<Face> Mesh::GatherInteriorFaces()
Kirill Terekhov's avatar
Kirill Terekhov committed
1260
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
1261
		ElementArray<Face> ret(this);
Kirill Terekhov's avatar
Kirill Terekhov committed
1262
1263
		if( GetMeshState() == Mesh::Serial )
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1264
1265
1266
1267
1268
			for(integer f = 0; f < FaceLastLocalID(); f++) if( isValidElement(FACE,f) )
			{
				Face ef = Face(this,ComposeHandle(FACE,f));
				if( ef->nbAdjElements(CELL) == 2 ) ret.push_back(ef);
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
1269
1270
1271
		}
		else
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1272
			for(integer f = 0; f < FaceLastLocalID(); f++) if( isValidElement(FACE,f) )
Kirill Terekhov's avatar
Kirill Terekhov committed
1273
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1274
1275
				Face ef = Face(this,ComposeHandle(FACE,f));
				ElementArray<Cell> cells = ef->getCells();
Kirill Terekhov's avatar
Kirill Terekhov committed
1276
				if( cells.size() == 2 && cells[0].GetStatus() != Element::Ghost && cells[1].GetStatus() != Element::Ghost)
Kirill Terekhov's avatar
Kirill Terekhov committed
1277
					ret.push_back(ef);
Kirill Terekhov's avatar
Kirill Terekhov committed
1278
1279
1280
1281
1282
1283
1284
			}
		}
		return ret;
	}

	Storage::integer Mesh::CountBoundaryFaces()
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
1285
1286
1287
1288
1289
1290
		integer ret = 0;
		for(integer f = 0; f < FaceLastLocalID(); f++) if( isValidElement(FACE,f) )
		{
			Face ef = Face(this,ComposeHandle(FACE,f));
			if( ef->Boundary() ) ret++;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
1291
1292
1293
1294
1295
		return ret;
	}

	Storage::integer Mesh::CountInteriorFaces()
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
1296
		integer ret = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1297
1298
		if( GetMeshState() == Mesh::Serial )
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1299
1300
1301
1302
1303
			for(integer f = 0; f < FaceLastLocalID(); f++) if( isValidElement(FACE,f) )
			{
				Face ef = Face(this,ComposeHandle(FACE,f));
				if( ef->nbAdjElements(CELL) == 2 ) ret++;
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
1304
1305
1306
		}
		else
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
1307
			for(integer f = 0; f < FaceLastLocalID(); f++) if( isValidElement(FACE,f) )
Kirill Terekhov's avatar
Kirill Terekhov committed
1308
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1309
1310
				Face ef = Face(this,ComposeHandle(FACE,f));
				ElementArray<Cell> cells = ef->getCells();
Kirill Terekhov's avatar
Kirill Terekhov committed
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
				if( cells.size() == 2 && cells[0].GetStatus() != Element::Ghost && cells[1].GetStatus() != Element::Ghost)
					ret++;
			}
		}

		return ret;
	}
	
	

	
Kirill Terekhov's avatar
Kirill Terekhov committed
1322
	void Element::CastRay(real * pos, real * dir, dynarray< std::pair<Element, real> , 16 > & hits) const
Kirill Terekhov's avatar
Kirill Terekhov committed
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
	{
		Mesh * mm = GetMeshLink();
		unsigned dim = mm->GetDimensions();
		Storage::real eps = mm->GetEpsilon();
		switch(GetElementType())
		{
			case NODE:
			{
				Storage::real_array coords = getAsNode()->Coords();
				Storage::real vec[3],ndir[3], lvec, ldir;
				for(unsigned k = 0; k < dim; k++)
				{
					vec[k] = (coords[k]-pos[k]);
					ndir[k] = dir[k];
				}
				lvec = vec_normalize(vec,dim);
				ldir = vec_normalize(ndir,dim);
				if( vec_dot_product(vec,ndir,dim) >= 1.0 - eps )
Kirill Terekhov's avatar
Kirill Terekhov committed
1341
					hits.push_back(std::make_pair(Element(mm,GetHandle()),lvec/ldir));
Kirill Terekhov's avatar
Kirill Terekhov committed
1342
1343
1344
1345
1346
1347
1348
1349
1350
			}	 
			break;
			case EDGE:
			{
				throw NotImplemented;