SubPolytope.cc 15 KB
Newer Older
1
2
3
4
5
#include "ElInfo.h"
#include "FixVec.h"
#include "SubElInfo.h"
#include "SubPolytope.h"

6
7
8
9
10
11
12
13
14
15
16
17
18
19
namespace AMDiS {

  bool SubPolytope::checkIntPoints()
  {
    ////////////////////////////////////////////////////////////////////////////
    //
    //  Do the points in intPoints lie on edges ? And are they inner points,
    //  i.e. they aren't vertices ?
    //  
    //  Return value: true  -  yes
    //                false -  no
    ////////////////////////////////////////////////////////////////////////////

    for (int i = 0; i < numIntPoints; i++) {
20
21
      int zeroCounter = 0;
      for (int j = 0; j < dim + 1; j++) {
22
23
24
	if (fabs((*intPoints)[i][j]) <= 1.e-15 ) { 
	  zeroCounter++; 
	}
25
      }
26

27
28
29
30
31
32
33
34
      /**
       *  Dimension 1
       *
       *  "Inner" points on edges aren't equal to 0.0 in any component.
       */
      if (dim == 1  &&  zeroCounter != 0 ) { 
	return false;
      }
35

36
37
38
39
40
41
42
43
      /**
       *  Dimension 2 
       *
       *  "Inner" points on edges are equal to 0.0 in exactly 1 component.
       */
      if (dim == 2  &&  zeroCounter != 1) {
	return false;
      }
44

45
46
47
48
49
50
51
52
      /**
       *  Dimension 3 
       *
       *  "Inner" points on edges are equal to 0.0 in exactly 2 components.
       */
      if (dim == 3  &&  zeroCounter != 2) { 
	return false;
      }
53
    }
54

55
56
    return true;
  }
57

58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
  SubPolytope::SubPolytope(const ElInfo *elInfo_, 
			   VectorOfFixVecs<DimVec<double> > *intPoints_,
			   int numIntPoints_,
			   const int &indexElVertInPol_)
    : elInfo(elInfo_),
      intPoints(intPoints_),
      numIntPoints(numIntPoints_)
  {
    FUNCNAME("SubPolytope::SubPolytope()");

    dim = (*intPoints_)[0].getSize() - 1;

    TEST_EXIT((dim == 1 && numIntPoints == 1) ||
	      (dim == 2 && numIntPoints == 2) ||
	      (dim == 3 && numIntPoints == 3) || 
	      (dim == 3 && numIntPoints == 4))
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
      ("invalid number of intersection points\n");

    /**
     *  Check whether the points in intPoints are really intersection points,
     *  i.e. lie on edges.
     */
    TEST_EXIT(checkIntPoints() == true)
      ("invalid intersection points - do not lie on edges\n");
    
    /*
     *  Create the subelements the polytope consists of.
     */
    switch ( dim ) {
    case 1:
      createSubElementPolytopeIsSubElement1D(indexElVertInPol_);
      break;
    case 2:
      createSubElementPolytopeIsSubElement2D3D();
      break;
    case 3: 
      if (numIntPoints == 3) {
	createSubElementPolytopeIsSubElement2D3D();
96
      } else {
97
98
99
100
101
102
103
	createSubElementsForPolytope3D(indexElVertInPol_);
      }
      break;
    default:
      ERROR_EXIT("invalid dimension\n");
      break;
    }
104
  }
105

106

107
108
109
110
111
112
113
  void SubPolytope::createSubElementPolytopeIsSubElement1D(int indexElVertInPol)
  {
    //**************************************************************************
    // The intersection of the one-dimensional element (interval) divides 
    // element into two subelements. indexElVertInPol indicates which 
    // subelement to take.
    //**************************************************************************
114

115
    FUNCNAME("SubPolytope::createSubElementPolytopeIsSubElement1D()");
116

117
    TEST_EXIT(dim == 1 && numIntPoints == 1)("invalid call of this routine\n");
118

119
    VectorOfFixVecs<DimVec<double> > *subElVertices = 
Thomas Witkowski's avatar
Thomas Witkowski committed
120
      new VectorOfFixVecs<DimVec<double> >(dim, dim + 1, NO_INIT);
121
    DimVec<double> vertex(dim, DEFAULT_VALUE, 1.0);
122
123

    /**
124
125
     *  Get the vertex which - with the intersection point in intPoints - forms
     *  a subelement of element.
126
     */
127
128
129
130
131
132
    if (indexElVertInPol == 0) {
      /**
       *  The vertex in element with barycentric coordinates (1,0) is in 
       *  subelement.
       */
      vertex[1] = 0.0;
133
    } else {
134
135
136
137
138
139
140
      /**
       *  The vertex in element with barycentric coordinates (0,1) is in 
       *  subelement.
       */
      vertex[0] = 0.0;
    }

141
    /**
142
143
144
145
146
147
     *  Collect the vertices of the subelement in subElVertices.
     *
     *  Note: The routines CompositeFEMOperator::getElementMatrix and 
     *        CompositeFEMOperator::getElementVector expect the first vertice 
     *        of a subelement to be a vertice of the corresponding element and
     *        not to be an intersection point.
148
     */
149
150
    (*subElVertices)[0] = vertex;
    (*subElVertices)[dim] = (*intPoints)[0];
151
152


153
154
155
    /**
     *  Create a new ElInfo for the subelement.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
156
    subElements.push_back( new SubElInfo(subElVertices, elInfo) );
157

158
159
    TEST_EXIT( subElements.size() == 1 )("error in creating subelements");
    numSubElements = 1;
160

Thomas Witkowski's avatar
Thomas Witkowski committed
161
    delete subElVertices;
162
  }
163
164


165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
  void SubPolytope::createSubElementPolytopeIsSubElement2D3D()
  {
    //**************************************************************************
    // The intersection with element produced dim intersection points.
    // Thus there is exactly one vertice in element which - with the 
    // intersection points - forms a subelement of element.
    // This routine determines this vertex and creates a new object of SubElInfo
    // for the subelement stored in subElements.
    //
    // How does it work ?
    // All intersection points lie on edges of element and aren't vertices. The
    // missing vertex is the intersection point of these edges.  
    // Lying on edges means, that each intersection point has exactly one 
    // component equal to 0.0. The missing vertex is equal to 0.0 in all these
    // components.
    //**************************************************************************

    FUNCNAME("SubPolytope::createSubElementPolytopeIsSubElement2D3D()");

    TEST_EXIT((dim == 2 && numIntPoints == 2) || 
	      (dim == 3 && numIntPoints == 3))
      ("invalid call of this routine\n");

    VectorOfFixVecs<DimVec<double> >*subElVertices = 
Thomas Witkowski's avatar
Thomas Witkowski committed
189
      new VectorOfFixVecs<DimVec<double> >(dim, dim + 1, NO_INIT);
190
    DimVec<double> vertex(dim, DEFAULT_VALUE, 1.0);
191

192
193
194
195
196
    /**
     *  Get the vertex which - with the intersection points intPoints - forms
     *  a subelement of element.
     */
    for (int i = 0; i < numIntPoints; i++) {
197
      for (int j = 0; j < dim + 1; j++) {
198
199
200
201
	if ( fabs((*intPoints)[i][j]) <= 1.e-15 ) {
	  vertex[j] = 0.0; 
	};
      }
202
203
    }

204
205
206
207
208
209
210
211
212
213
214
215
    /**
     *  Collect the vertices of the subelement in subElVertices.
     *
     *  Note: The routines CompositeFEMOperator::getElementMatrix and 
     *        CompositeFEMOperator::getElementVector expect the first vertice 
     *        of a subelement to be a vertice of the corresponding element and
     *        not to be an intersection point.
     */
    (*subElVertices)[0] = vertex;
    for (int i = 0; i < numIntPoints; i++) {
      (*subElVertices)[i+1] = (*intPoints)[i];
    }
216
217


218
219
220
    /**
     *  Create a new ElInfo for the subelement.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
221
    subElements.push_back( new SubElInfo(subElVertices, elInfo) );
222

223
224
    TEST_EXIT( subElements.size() == 1 )("error in creating subelements");
    numSubElements = 1;
225

Thomas Witkowski's avatar
Thomas Witkowski committed
226
    delete subElVertices;
227
  }
228
229


230
231
  int SubPolytope::getIndexSecondFaceIntPoint0(int indexFirstFace, int dim)
  {
232
    for (int i = 0; i < dim + 1; i++) {
233
234
235
      if ( fabs((*intPoints)[0][i]) <= 1.e-15  &&  i != indexFirstFace ) {
	return i;
      }
236
    }
237

238
239
240
    ERROR_EXIT("couldn't determine the second face for IntPoint0\n");
    return -1;
  }
241

242
 
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
  void SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
  {
    //**************************************************************************
    // The intersection with element produced four intersection points. Thus the
    // intersection doesn't induce a subelement. This routine divides the
    // subpolytope given by the intersection into three subelements.
    //
    // How does it work ?
    // The intersection points and two vertices of element form a subpolytope
    // of element. First of all, we determine these vertices, and call them
    // A and B. Then we sort the intersection points (S_0, S_1, S_2, S_3)
    // in the following way:
    // S_0 is the first intersection point in the creation-list /ref intPoints_.
    // A and S_0 lie in the face of element opposite to B.
    // S_0 and S_1 are neighbours of A.
    // S_2 is opposite to S_0 in the intersection plane.
    // S_1 and S_3 are neighbours of A.
    // Then the subelements are:
    // A - B - S_0 - S_1 ,  B - S_0 - S_1 - S_2 ,  B - S_0 - S_2 - S_3 .
    //
    // The index of one vertex of element that is in subpolytope is handed to
    // this routine. The subpolytope is determined uniquely by this vertex and
    // the intersection points.
    //**************************************************************************

    FUNCNAME("SubPolytope::createSubElementForPolytope3D");

    TEST_EXIT(dim == 3 && numIntPoints == 4)("invalid call of this routine\n");

    TEST_EXIT(0 <= indexElVertInPol1  &&  indexElVertInPol1 <= 3)
      ("invalid index for vertex of a tetrahedron");

    VectorOfFixVecs<DimVec<double> > *subElVertices = 
Thomas Witkowski's avatar
Thomas Witkowski committed
276
      new VectorOfFixVecs<DimVec<double> >(dim, dim + 1, NO_INIT);
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
    DimVec<double> vertexA(dim, DEFAULT_VALUE, 0.0);
    DimVec<double> vertexB(dim, DEFAULT_VALUE, 0.0);

    int indexElVertInPol2 = 0;  // index of second vertex of element lying in 
    // subpolytope
    // The vertices of element (3D -> tetrahedron) are
    // indexed as usual:
    // 0: (1,0,0,0), 1: (0,1,0,0), 2: (0,0,1,0), 
    // 3: (0,0,0,1)
    bool intPointOnEdge[4][4] = {{false, false, false, false},
				 {false, false, false, false},
				 {false, false, false, false},
				 {false, false, false, false}};
    // /ref intPointOnEdge[i][j] indicates whether there
    // is an intersection point on the edge from  
    // vertice i to vertice j :
    // false : no intersection point
    // true: there is an intersection point
    int indexEdge[2];       // For a vertex lying on an edge indexEdge 
    // stores the indices of the two barycentric 
    // coordinates which are not equal to zero.
    int indexA = 0;            
    int indexB = 0;            
    int indexS_0 = 0;
    int indexS_1 = 0;
    int indexS_2 = 0;
    int indexS_3 = 0;
    int indexSecondFaceIntPoint0 = 0;

    /**
     *  Get the second vertex of element lying in subpolytope.
     *
     *  There is exactly one vertex of element which - with vertex 
     *  indexElVertInPol1 and the intersection points intPoints -
     *  forms a subpolytope of element. It is the vertex adjacent with
     *  indexElVertInPol1 whose common edge with indexElVertInPol1
     *  doesn't contain an intersection point.
     */

    // Get the edges including the intersection points.
    for (int i = 0; i < numIntPoints; i++) {
      int k = 0;
319
      for (int j = 0; j < dim + 1; j++) {
320
321
322
323
	if (fabs((*intPoints)[i][j]) > 1.e-15 ) {
	  indexEdge[k] = j;
	  k++;
	}
324
      }
325
326
      intPointOnEdge[indexEdge[0]][indexEdge[1]] = true;
      intPointOnEdge[indexEdge[1]][indexEdge[0]] = true;
327
328
    }

329
330
331
    // Get the vertex of element adjacent with indexElVertInPol1 whose
    // common edge with indexElVertInPol1 doesn't contain an
    // intersection point, and store it in indexElVertInPol2.
332
    for (int i = 0; i < dim + 1; i++) {
333
334
335
336
337
      if (intPointOnEdge[indexElVertInPol1][i] == false  &&  
	  i != indexElVertInPol1 ) {
	indexElVertInPol2 = i;
	break;
      }
338
339
    }

340
341
342
343
344
345
346
    /**
     *  Determine A and B, so that intPoint0 is a neighbour of A. 
     *
     *  In the subpolytope A and two intersection points lie in the face
     *  opposite to B. And B and the other two intersection points lie in the
     *  face opposite to A.
     */
347

348
    if (fabs((*intPoints)[0][indexElVertInPol1]) <= 1.e-15) {
349

350
351
      // (*intPoints)[0] lies in the face opposite to vertex 
      // /ref indexElVertInPol1.
352

353
354
355
      indexA = indexElVertInPol2;
      indexB = indexElVertInPol1;
    } else if (fabs((*intPoints)[0][indexElVertInPol2]) <= 1.e-15) {
356

357
358
      // (*intPoints)[0] lies in the face opposite to vertex
      // /ref indexElVertInPol2.
359

360
361
362
363
364
      indexA = indexElVertInPol1;
      indexB = indexElVertInPol2;
    } else {
      ERROR_EXIT("couldn't determine A and B\n");
    }
365

366
367
368
    /**
     *  Sort the intersection points.
     */
369

370
371
    // (*intPoints)[0] is a neighbour of A (A has been constructed this way).
    indexS_0 = 0;
372

373
    if (fabs((*intPoints)[1][indexB]) <= 1.e-15) {
374

375
376
      // (*intPoints)[1] lies in the face opposite to B, thus is a neighbour
      // of A.
377

378
      indexS_1 = 1;
379

380
      indexSecondFaceIntPoint0 = getIndexSecondFaceIntPoint0(indexB, dim);
381

382
      if (fabs((*intPoints)[2][indexSecondFaceIntPoint0]) <= 1.e-15) {
383

384
	// (*intPoints)[2] is neighbour of (*intPoints)[0] 
385

386
387
388
	indexS_2 = 3;
	indexS_3 = 2;
      } else {
389

390
391
	// (*intPoints)[2] is opposite to (*intPoints)[0] in the intersection
	// plane
392

393
394
395
396
	indexS_2 = 2;
	indexS_3 = 3;
      }
    } else if (fabs((*intPoints)[1][indexA]) <= 1.e-15) {
397

398
      // (*intPoints)[1] lies in the face opposite to A
399

400
      indexSecondFaceIntPoint0 = getIndexSecondFaceIntPoint0(indexB, dim);
401

402
      if (fabs((*intPoints)[1][indexSecondFaceIntPoint0]) <= 1.e-15) {
403

404
405
	// (*intPoints)[1] is neighbour of (*intPoints)[0], but isn't 
	// neighbour of A
406

407
	indexS_3 = 1;
408

409
	if (fabs((*intPoints)[2][indexB]) <= 1.e-15) {
410

411
	  // (*intPoints)[2] is neighbour of (*intPoints)[0] and neighbour of A
412

413
414
415
	  indexS_1 = 2;
	  indexS_2 = 3;
	} else {
416

417
418
	  // (*intPoints)[2] is opposite to (*intPoints)[0] in the intersection
	  // plane
419

420
421
422
423
	  indexS_1 = 3;
	  indexS_2 = 2;
	}
      } else {
424

425
426
	// (*intPoints)[1] isn't neighbour of (*intPoints)[0], thus lies opposite 
	// to (*intPoints)[0] in the intersection plane
427

428
	indexS_2 = 1;
429

430
	if (fabs((*intPoints)[2][indexB]) <= 1.e-15) {
431

432
	  // (*intPoints)[2] is neighbour of A
433

434
435
436
	  indexS_1 = 2;
	  indexS_3 = 3;
	} else {
437

438
	  // (*intPoints)[2] isn't neighbour of A
439

440
441
442
	  indexS_1 = 3;
	  indexS_3 = 2;
	}
443
      }
444
445
    } else {
      ERROR_EXIT("IntPoint1 isn't either part of the face opposite to A nor part of the face opposite to B\n");
446
    }
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467

    /**
     *  For each subelement: Collect the vertices of the subelement in 
     *  subElVertices, create a new SubElInfo for the subelement and
     *  store it in subElements.
     *
     *  Note: The routines CompositeFEMOperator::getElementMatrix and 
     *        CompositeFEMOperator::getElementVector expect the first vertice 
     *        of a subelement to be a vertice of the corresponding element and
     *        not to be an intersection point.
     */

    // Create vertex A and vertex B.
    vertexA[indexA] = 1.0;
    vertexB[indexB] = 1.0;

    // Subelement 1: A - B - S_0 - S_1
    (*subElVertices)[0] = vertexA;
    (*subElVertices)[1] = vertexB;
    (*subElVertices)[2] = (*intPoints)[indexS_0];
    (*subElVertices)[3] = (*intPoints)[indexS_1];
Thomas Witkowski's avatar
Thomas Witkowski committed
468
    subElements.push_back( new SubElInfo(subElVertices, elInfo) );
469
470
471
472
473
474

    // Subelement 2: B - S_0 - S_1 - S_2
    (*subElVertices)[0] = vertexB;
    (*subElVertices)[1] = (*intPoints)[indexS_0];
    (*subElVertices)[2] = (*intPoints)[indexS_1];
    (*subElVertices)[3] = (*intPoints)[indexS_2];
Thomas Witkowski's avatar
Thomas Witkowski committed
475
    subElements.push_back( new SubElInfo(subElVertices, elInfo) );
476
477
478
479
480
481

    // Subelement 3: B - S_0 - S_2 - S_3
    (*subElVertices)[0] = vertexB;
    (*subElVertices)[1] = (*intPoints)[indexS_0];
    (*subElVertices)[2] = (*intPoints)[indexS_2];
    (*subElVertices)[3] = (*intPoints)[indexS_3];
Thomas Witkowski's avatar
Thomas Witkowski committed
482
    subElements.push_back( new SubElInfo(subElVertices, elInfo) );
483
484
485
486

    TEST_EXIT( subElements.size() == 3 )("error in creating subelements");
    numSubElements = 3;

Thomas Witkowski's avatar
Thomas Witkowski committed
487
    delete subElVertices;
488
489
490
  }

}