MFEM  v3.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
geom.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.googlecode.com.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "fem.hpp"
13 
14 namespace mfem
15 {
16 
17 const char *Geometry::Name[NumGeom] =
18 { "Point", "Segment", "Triangle", "Square", "Tetrahedron", "Cube" };
19 
20 const double Geometry::Volume[NumGeom] =
21 { 1.0, 1.0, 0.5, 1.0, 1./6, 1.0 };
22 
23 Geometry::Geometry()
24 {
25  // Vertices for Geometry::POINT
26  GeomVert[0] = NULL; // No vertices, dimension is 0
27 
28  // Vertices for Geometry::SEGMENT
29  GeomVert[1] = new IntegrationRule(2);
30 
31  GeomVert[1]->IntPoint(0).x = 0.0;
32  GeomVert[1]->IntPoint(1).x = 1.0;
33 
34  // Vertices for Geometry::TRIANGLE
35  GeomVert[2] = new IntegrationRule(3);
36 
37  GeomVert[2]->IntPoint(0).x = 0.0;
38  GeomVert[2]->IntPoint(0).y = 0.0;
39 
40  GeomVert[2]->IntPoint(1).x = 1.0;
41  GeomVert[2]->IntPoint(1).y = 0.0;
42 
43  GeomVert[2]->IntPoint(2).x = 0.0;
44  GeomVert[2]->IntPoint(2).y = 1.0;
45 
46  // Vertices for Geometry::SQUARE
47  GeomVert[3] = new IntegrationRule(4);
48 
49  GeomVert[3]->IntPoint(0).x = 0.0;
50  GeomVert[3]->IntPoint(0).y = 0.0;
51 
52  GeomVert[3]->IntPoint(1).x = 1.0;
53  GeomVert[3]->IntPoint(1).y = 0.0;
54 
55  GeomVert[3]->IntPoint(2).x = 1.0;
56  GeomVert[3]->IntPoint(2).y = 1.0;
57 
58  GeomVert[3]->IntPoint(3).x = 0.0;
59  GeomVert[3]->IntPoint(3).y = 1.0;
60 
61  // Vertices for Geometry::TETRAHEDRON
62  GeomVert[4] = new IntegrationRule(4);
63  GeomVert[4]->IntPoint(0).x = 0.0;
64  GeomVert[4]->IntPoint(0).y = 0.0;
65  GeomVert[4]->IntPoint(0).z = 0.0;
66 
67  GeomVert[4]->IntPoint(1).x = 1.0;
68  GeomVert[4]->IntPoint(1).y = 0.0;
69  GeomVert[4]->IntPoint(1).z = 0.0;
70 
71  GeomVert[4]->IntPoint(2).x = 0.0;
72  GeomVert[4]->IntPoint(2).y = 1.0;
73  GeomVert[4]->IntPoint(2).z = 0.0;
74 
75  GeomVert[4]->IntPoint(3).x = 0.0;
76  GeomVert[4]->IntPoint(3).y = 0.0;
77  GeomVert[4]->IntPoint(3).z = 1.0;
78 
79  // Vertices for Geometry::CUBE
80  GeomVert[5] = new IntegrationRule(8);
81 
82  GeomVert[5]->IntPoint(0).x = 0.0;
83  GeomVert[5]->IntPoint(0).y = 0.0;
84  GeomVert[5]->IntPoint(0).z = 0.0;
85 
86  GeomVert[5]->IntPoint(1).x = 1.0;
87  GeomVert[5]->IntPoint(1).y = 0.0;
88  GeomVert[5]->IntPoint(1).z = 0.0;
89 
90  GeomVert[5]->IntPoint(2).x = 1.0;
91  GeomVert[5]->IntPoint(2).y = 1.0;
92  GeomVert[5]->IntPoint(2).z = 0.0;
93 
94  GeomVert[5]->IntPoint(3).x = 0.0;
95  GeomVert[5]->IntPoint(3).y = 1.0;
96  GeomVert[5]->IntPoint(3).z = 0.0;
97 
98  GeomVert[5]->IntPoint(4).x = 0.0;
99  GeomVert[5]->IntPoint(4).y = 0.0;
100  GeomVert[5]->IntPoint(4).z = 1.0;
101 
102  GeomVert[5]->IntPoint(5).x = 1.0;
103  GeomVert[5]->IntPoint(5).y = 0.0;
104  GeomVert[5]->IntPoint(5).z = 1.0;
105 
106  GeomVert[5]->IntPoint(6).x = 1.0;
107  GeomVert[5]->IntPoint(6).y = 1.0;
108  GeomVert[5]->IntPoint(6).z = 1.0;
109 
110  GeomVert[5]->IntPoint(7).x = 0.0;
111  GeomVert[5]->IntPoint(7).y = 1.0;
112  GeomVert[5]->IntPoint(7).z = 1.0;
113 
114  GeomCenter[POINT].x = 0.0;
115  GeomCenter[POINT].y = 0.0;
116  GeomCenter[POINT].z = 0.0;
117 
118  GeomCenter[SEGMENT].x = 0.5;
119  GeomCenter[SEGMENT].y = 0.0;
120  GeomCenter[SEGMENT].z = 0.0;
121 
122  GeomCenter[TRIANGLE].x = 1.0 / 3.0;
123  GeomCenter[TRIANGLE].y = 1.0 / 3.0;
124  GeomCenter[TRIANGLE].z = 0.0;
125 
126  GeomCenter[SQUARE].x = 0.5;
127  GeomCenter[SQUARE].y = 0.5;
128  GeomCenter[SQUARE].z = 0.0;
129 
130  GeomCenter[TETRAHEDRON].x = 0.25;
131  GeomCenter[TETRAHEDRON].y = 0.25;
132  GeomCenter[TETRAHEDRON].z = 0.25;
133 
134  GeomCenter[CUBE].x = 0.5;
135  GeomCenter[CUBE].y = 0.5;
136  GeomCenter[CUBE].z = 0.5;
137 
138  PerfGeomToGeomJac[POINT] = NULL;
139  PerfGeomToGeomJac[SEGMENT] = NULL;
140  PerfGeomToGeomJac[TRIANGLE] = new DenseMatrix(2);
141  PerfGeomToGeomJac[SQUARE] = NULL;
142  PerfGeomToGeomJac[TETRAHEDRON] = new DenseMatrix(3);
143  PerfGeomToGeomJac[CUBE] = NULL;
144 
145  {
146  Linear2DFiniteElement TriFE;
147  IsoparametricTransformation tri_T;
148  tri_T.SetFE(&TriFE);
149  GetPerfPointMat (TRIANGLE, tri_T.GetPointMat());
150  tri_T.SetIntPoint(&GeomCenter[TRIANGLE]);
151  CalcInverse(tri_T.Jacobian(), *PerfGeomToGeomJac[TRIANGLE]);
152  }
153  {
154  Linear3DFiniteElement TetFE;
155  IsoparametricTransformation tet_T;
156  tet_T.SetFE(&TetFE);
157  GetPerfPointMat (TETRAHEDRON, tet_T.GetPointMat());
158  tet_T.SetIntPoint(&GeomCenter[TETRAHEDRON]);
159  CalcInverse(tet_T.Jacobian(), *PerfGeomToGeomJac[TETRAHEDRON]);
160  }
161 }
162 
163 Geometry::~Geometry()
164 {
165  for (int i = 0; i < NumGeom; i++)
166  {
167  delete PerfGeomToGeomJac[i];
168  delete GeomVert[i];
169  }
170 }
171 
172 const IntegrationRule * Geometry::GetVertices(int GeomType)
173 {
174  switch (GeomType)
175  {
176  case Geometry::POINT: return GeomVert[0];
177  case Geometry::SEGMENT: return GeomVert[1];
178  case Geometry::TRIANGLE: return GeomVert[2];
179  case Geometry::SQUARE: return GeomVert[3];
180  case Geometry::TETRAHEDRON: return GeomVert[4];
181  case Geometry::CUBE: return GeomVert[5];
182  default:
183  mfem_error ("Geometry::GetVertices(...)");
184  }
185  // make some compilers happy.
186  return GeomVert[0];
187 }
188 
189 void Geometry::GetPerfPointMat(int GeomType, DenseMatrix &pm)
190 {
191  switch (GeomType)
192  {
193  case Geometry::SEGMENT:
194  {
195  pm.SetSize (1, 2);
196  pm(0,0) = 0.0;
197  pm(0,1) = 1.0;
198  }
199  break;
200 
201  case Geometry::TRIANGLE:
202  {
203  pm.SetSize (2, 3);
204  pm(0,0) = 0.0; pm(1,0) = 0.0;
205  pm(0,1) = 1.0; pm(1,1) = 0.0;
206  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676;
207  }
208  break;
209 
210  case Geometry::SQUARE:
211  {
212  pm.SetSize (2, 4);
213  pm(0,0) = 0.0; pm(1,0) = 0.0;
214  pm(0,1) = 1.0; pm(1,1) = 0.0;
215  pm(0,2) = 1.0; pm(1,2) = 1.0;
216  pm(0,3) = 0.0; pm(1,3) = 1.0;
217  }
218  break;
219 
220  case Geometry::TETRAHEDRON:
221  {
222  pm.SetSize (3, 4);
223  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
224  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
225  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676; pm(2,2) = 0.0;
226  pm(0,3) = 0.5; pm(1,3) = 0.28867513459481288225;
227  pm(2,3) = 0.81649658092772603273;
228  }
229  break;
230 
231  case Geometry::CUBE:
232  {
233  pm.SetSize (3, 8);
234  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
235  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
236  pm(0,2) = 1.0; pm(1,2) = 1.0; pm(2,2) = 0.0;
237  pm(0,3) = 0.0; pm(1,3) = 1.0; pm(2,3) = 0.0;
238  pm(0,4) = 0.0; pm(1,4) = 0.0; pm(2,4) = 1.0;
239  pm(0,5) = 1.0; pm(1,5) = 0.0; pm(2,5) = 1.0;
240  pm(0,6) = 1.0; pm(1,6) = 1.0; pm(2,6) = 1.0;
241  pm(0,7) = 0.0; pm(1,7) = 1.0; pm(2,7) = 1.0;
242  }
243  break;
244 
245  default:
246  mfem_error ("Geometry::GetPerfPointMat (...)");
247  }
248 }
249 
250 void Geometry::JacToPerfJac(int GeomType, const DenseMatrix &J,
251  DenseMatrix &PJ) const
252 {
253  if (PerfGeomToGeomJac[GeomType])
254  {
255  Mult(J, *PerfGeomToGeomJac[GeomType], PJ);
256  }
257  else
258  {
259  PJ = J;
260  }
261 }
262 
263 const int Geometry::NumGeom;
264 
265 const int Geometry::NumBdrArray[] = { 0, 2, 3, 4, 4, 6 };
266 
267 Geometry Geometries;
268 
269 
270 GeometryRefiner::GeometryRefiner()
271 {
272  type = 0;
273  for (int i = 0; i < Geometry::NumGeom; i++)
274  {
275  RGeom[i] = NULL;
276  IntPts[i] = NULL;
277  }
278 }
279 
280 GeometryRefiner::~GeometryRefiner()
281 {
282  for (int i = 0; i < Geometry::NumGeom; i++)
283  {
284  delete RGeom[i];
285  delete IntPts[i];
286  }
287 }
288 
289 RefinedGeometry * GeometryRefiner::Refine (int Geom, int Times, int ETimes)
290 {
291  int i, j, k, l;
292 
293  const double *cp = NULL;
294  if (type)
295  cp = poly1d.ClosedPoints(Times);
296 
297  switch (Geom)
298  {
299  case Geometry::SEGMENT:
300  {
301  const int g = Geometry::SEGMENT;
302  if (RGeom[g] != NULL && RGeom[g]->Times == Times)
303  return RGeom[g];
304  delete RGeom[g];
305  RGeom[g] = new RefinedGeometry(Times+1, 2*Times, 0);
306  RGeom[g]->Times = Times;
307  RGeom[g]->ETimes = 0;
308  for (i = 0; i <= Times; i++)
309  {
310  IntegrationPoint &ip = RGeom[g]->RefPts.IntPoint(i);
311  ip.x = (type == 0) ? double(i) / Times : cp[i];
312  }
313  Array<int> &G = RGeom[g]->RefGeoms;
314  for (i = 0; i < Times; i++)
315  {
316  G[2*i+0] = i;
317  G[2*i+1] = i+1;
318  }
319 
320  return RGeom[g];
321  }
322 
323  case Geometry::TRIANGLE:
324  {
325  if (RGeom[2] != NULL && RGeom[2]->Times == Times &&
326  RGeom[2]->ETimes == ETimes)
327  return RGeom[2];
328 
329  if (RGeom[2] != NULL)
330  delete RGeom[2];
331  RGeom[2] = new RefinedGeometry ((Times+1)*(Times+2)/2, 3*Times*Times,
332  3*Times*(ETimes+1));
333  RGeom[2]->Times = Times;
334  RGeom[2]->ETimes = ETimes;
335  for (k = j = 0; j <= Times; j++)
336  for (i = 0; i <= Times-j; i++, k++)
337  {
338  IntegrationPoint &ip = RGeom[2]->RefPts.IntPoint(k);
339  if (type == 0)
340  {
341  ip.x = double(i) / Times;
342  ip.y = double(j) / Times;
343  }
344  else
345  {
346  ip.x = cp[i]/(cp[i] + cp[j] + cp[Times-i-j]);
347  ip.y = cp[j]/(cp[i] + cp[j] + cp[Times-i-j]);
348  }
349  }
350  Array<int> &G = RGeom[2]->RefGeoms;
351  for (l = k = j = 0; j < Times; j++, k++)
352  for (i = 0; i < Times-j; i++, k++)
353  {
354  G[l++] = k;
355  G[l++] = k+1;
356  G[l++] = k+Times-j+1;
357  if (i+j+1 < Times)
358  {
359  G[l++] = k+1;
360  G[l++] = k+Times-j+2;
361  G[l++] = k+Times-j+1;
362  }
363  }
364  Array<int> &E = RGeom[2]->RefEdges;
365  // horizontal edges
366  for (l = k = 0; k < Times; k += Times/ETimes)
367  {
368  j = k*(Times+1)-((k-1)*k)/2;
369  for (i = 0; i < Times-k; i++)
370  {
371  E[l++] = j; j++;
372  E[l++] = j;
373  }
374  }
375  // diagonal edges
376  for (k = Times; k > 0; k -= Times/ETimes)
377  {
378  j = k;
379  for (i = 0; i < k; i++)
380  {
381  E[l++] = j; j += Times-i;
382  E[l++] = j;
383  }
384  }
385  // vertical edges
386  for (k = 0; k < Times; k += Times/ETimes)
387  {
388  j = k;
389  for (i = 0; i < Times-k; i++)
390  {
391  E[l++] = j; j += Times-i+1;
392  E[l++] = j;
393  }
394  }
395 
396  return RGeom[2];
397  }
398 
399  case Geometry::SQUARE:
400  {
401  if (RGeom[3] != NULL && RGeom[3]->Times == Times &&
402  RGeom[3]->ETimes == ETimes)
403  return RGeom[3];
404 
405  if (RGeom[3] != NULL)
406  delete RGeom[3];
407  RGeom[3] = new RefinedGeometry ((Times+1)*(Times+1), 4*Times*Times,
408  4*(ETimes+1)*Times);
409  RGeom[3]->Times = Times;
410  RGeom[3]->ETimes = ETimes;
411  for (k = j = 0; j <= Times; j++)
412  for (i = 0; i <= Times; i++, k++)
413  {
414  IntegrationPoint &ip = RGeom[3]->RefPts.IntPoint(k);
415  if (type == 0)
416  {
417  ip.x = double(i) / Times;
418  ip.y = double(j) / Times;
419  }
420  else
421  {
422  ip.x = cp[i];
423  ip.y = cp[j];
424  }
425  }
426  Array<int> &G = RGeom[3]->RefGeoms;
427  for (l = k = j = 0; j < Times; j++, k++)
428  for (i = 0; i < Times; i++, k++)
429  {
430  G[l++] = k;
431  G[l++] = k+1;
432  G[l++] = k+Times+2;
433  G[l++] = k+Times+1;
434  }
435  Array<int> &E = RGeom[3]->RefEdges;
436  // horizontal edges
437  for (l = k = 0; k <= Times; k += Times/ETimes)
438  {
439  for (i = 0, j = k*(Times+1); i < Times; i++)
440  {
441  E[l++] = j; j++;
442  E[l++] = j;
443  }
444  }
445  // vertical edges (in right-to-left order)
446  for (k = Times; k >= 0; k -= Times/ETimes)
447  {
448  for (i = 0, j = k; i < Times; i++)
449  {
450  E[l++] = j; j += Times+1;
451  E[l++] = j;
452  }
453  }
454 
455  return RGeom[3];
456  }
457 
458  case Geometry::CUBE:
459  {
460  const int g = Geometry::CUBE;
461  if (RGeom[g] != NULL && RGeom[g]->Times == Times &&
462  RGeom[g]->ETimes == ETimes)
463  return RGeom[g];
464 
465  if (RGeom[g] != NULL)
466  delete RGeom[g];
467  RGeom[g] = new RefinedGeometry ((Times+1)*(Times+1)*(Times+1),
468  8*Times*Times*Times, 0);
469  RGeom[g]->Times = Times;
470  RGeom[g]->ETimes = ETimes;
471  for (l = k = 0; k <= Times; k++)
472  for (j = 0; j <= Times; j++)
473  for (i = 0; i <= Times; i++, l++)
474  {
475  IntegrationPoint &ip = RGeom[g]->RefPts.IntPoint(l);
476  if (type == 0)
477  {
478  ip.x = double(i) / Times;
479  ip.y = double(j) / Times;
480  ip.z = double(k) / Times;
481  }
482  else
483  {
484  ip.x = cp[i];
485  ip.y = cp[j];
486  ip.z = cp[k];
487  }
488  }
489  Array<int> &G = RGeom[g]->RefGeoms;
490  for (l = k = 0; k < Times; k++)
491  for (j = 0; j < Times; j++)
492  for (i = 0; i < Times; i++)
493  {
494  G[l++] = i+0 + (j+0 + (k+0) * (Times+1)) * (Times+1);
495  G[l++] = i+1 + (j+0 + (k+0) * (Times+1)) * (Times+1);
496  G[l++] = i+1 + (j+1 + (k+0) * (Times+1)) * (Times+1);
497  G[l++] = i+0 + (j+1 + (k+0) * (Times+1)) * (Times+1);
498  G[l++] = i+0 + (j+0 + (k+1) * (Times+1)) * (Times+1);
499  G[l++] = i+1 + (j+0 + (k+1) * (Times+1)) * (Times+1);
500  G[l++] = i+1 + (j+1 + (k+1) * (Times+1)) * (Times+1);
501  G[l++] = i+0 + (j+1 + (k+1) * (Times+1)) * (Times+1);
502  }
503 
504  return RGeom[g];
505  }
506 
507  case Geometry::TETRAHEDRON:
508  {
509  const int g = Geometry::TETRAHEDRON;
510  if (RGeom[g] != NULL && RGeom[g]->Times == Times &&
511  RGeom[g]->ETimes == ETimes)
512  return RGeom[g];
513 
514  if (RGeom[g] != NULL)
515  delete RGeom[g];
516 
517  // subdivide the tetrahedron with vertices
518  // (0,0,0), (0,0,1), (1,1,1), (0,1,1)
519 
520  // vertices: 0 <= i <= j <= k <= Times
521  // (3-combination with repetitions)
522  // number of vertices: (n+3)*(n+2)*(n+1)/6, n = Times
523 
524  // elements: the vertices are: v1=(i,j,k), v2=v1+u1, v3=v2+u2, v4=v3+u3
525  // where 0 <= i <= j <= k <= n-1 and
526  // u1,u2,u3 is a permutation of (1,0,0),(0,1,0),(0,0,1)
527  // such that all v2,v3,v4 have non-decreasing components
528  // number of elements: n^3
529 
530  const int n = Times;
531  RGeom[g] = new RefinedGeometry((n+3)*(n+2)*(n+1)/6, 4*n*n*n, 0);
532  // enumerate and define the vertices
533  Array<int> vi((n+1)*(n+1)*(n+1));
534  vi = -1;
535  int m = 0;
536  for (k = 0; k <= n; k++)
537  for (j = 0; j <= k; j++)
538  for (i = 0; i <= j; i++)
539  {
540  IntegrationPoint &ip = RGeom[g]->RefPts.IntPoint(m);
541  // map the coordinates to the reference tetrahedron
542  // (0,0,0) -> (0,0,0)
543  // (0,0,1) -> (1,0,0)
544  // (1,1,1) -> (0,1,0)
545  // (0,1,1) -> (0,0,1)
546  if (type == 0)
547  {
548  ip.x = double(k - j) / n;
549  ip.y = double(i) / n;
550  ip.z = double(j - i) / n;
551  }
552  else
553  {
554  double w = cp[k-j] + cp[i] + cp[j-i] + cp[Times-k];
555  ip.x = cp[k-j]/w;
556  ip.y = cp[i]/w;
557  ip.z = cp[j-i]/w;
558  }
559  l = i + (j + k * (n+1)) * (n+1);
560  vi[l] = m;
561  m++;
562  }
563  if (m != (n+3)*(n+2)*(n+1)/6)
564  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #1");
565  // elements
566  Array<int> &G = RGeom[g]->RefGeoms;
567  m = 0;
568  for (k = 0; k < n; k++)
569  for (j = 0; j <= k; j++)
570  for (i = 0; i <= j; i++)
571  {
572  // the ordering of the vertices is chosen to ensure:
573  // 1) correct orientation
574  // 2) the x,y,z edges are in the set of edges
575  // {(0,1),(2,3), (0,2),(1,3)}
576  // (goal is to ensure that subsequent refinement using
577  // this procedure preserves the six tetrahedral shapes)
578 
579  // zyx: (i,j,k)-(i,j,k+1)-(i+1,j+1,k+1)-(i,j+1,k+1)
580  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
581  G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
582  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
583  G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
584  if (j < k)
585  {
586  // yzx: (i,j,k)-(i+1,j+1,k+1)-(i,j+1,k)-(i,j+1,k+1)
587  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
588  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
589  G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
590  G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
591  // yxz: (i,j,k)-(i,j+1,k)-(i+1,j+1,k+1)-(i+1,j+1,k)
592  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
593  G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
594  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
595  G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
596  }
597  if (i < j)
598  {
599  // xzy: (i,j,k)-(i+1,j,k)-(i+1,j+1,k+1)-(i+1,j,k+1)
600  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
601  G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
602  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
603  G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
604  if (j < k)
605  {
606  // xyz: (i,j,k)-(i+1,j+1,k+1)-(i+1,j,k)-(i+1,j+1,k)
607  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
608  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
609  G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
610  G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
611  }
612  // zxy: (i,j,k)-(i+1,j+1,k+1)-(i,j,k+1)-(i+1,j,k+1)
613  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
614  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
615  G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
616  G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
617  }
618  }
619  if (m != 4*n*n*n)
620  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #2");
621  for (i = 0; i < m; i++)
622  if (G[i] < 0)
623  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #3");
624 
625  return RGeom[g];
626  }
627 
628  default:
629 
630  return RGeom[0];
631  }
632 }
633 
634 const IntegrationRule *GeometryRefiner::RefineInterior(int Geom, int Times)
635 {
636  int g = Geom;
637 
638  switch (g)
639  {
640  case Geometry::SEGMENT:
641  {
642  if (Times < 2)
643  return NULL;
644  if (IntPts[g] == NULL || IntPts[g]->GetNPoints() != Times-1)
645  {
646  delete IntPts[g];
647  IntPts[g] = new IntegrationRule(Times-1);
648  for (int i = 1; i < Times; i++)
649  {
650  IntegrationPoint &ip = IntPts[g]->IntPoint(i-1);
651  ip.x = double(i) / Times;
652  ip.y = ip.z = 0.0;
653  }
654  }
655  }
656  break;
657 
658  case Geometry::TRIANGLE:
659  {
660  if (Times < 3)
661  return NULL;
662  if (IntPts[g] == NULL ||
663  IntPts[g]->GetNPoints() != ((Times-1)*(Times-2))/2)
664  {
665  delete IntPts[g];
666  IntPts[g] = new IntegrationRule(((Times-1)*(Times-2))/2);
667  for (int k = 0, j = 1; j < Times-1; j++)
668  for (int i = 1; i < Times-j; i++, k++)
669  {
670  IntegrationPoint &ip = IntPts[g]->IntPoint(k);
671  ip.x = double(i) / Times;
672  ip.y = double(j) / Times;
673  ip.z = 0.0;
674  }
675  }
676  }
677  break;
678 
679  case Geometry::SQUARE:
680  {
681  if (Times < 2)
682  return NULL;
683  if (IntPts[g] == NULL || IntPts[g]->GetNPoints() != (Times-1)*(Times-1))
684  {
685  delete IntPts[g];
686  IntPts[g] = new IntegrationRule((Times-1)*(Times-1));
687  for (int k = 0, j = 1; j < Times; j++)
688  for (int i = 1; i < Times; i++, k++)
689  {
690  IntegrationPoint &ip = IntPts[g]->IntPoint(k);
691  ip.x = double(i) / Times;
692  ip.y = double(j) / Times;
693  ip.z = 0.0;
694  }
695  }
696  }
697  break;
698 
699  default:
700  mfem_error("GeometryRefiner::RefineInterior(...)");
701  }
702 
703  return IntPts[g];
704 }
705 
706 GeometryRefiner GlobGeometryRefiner;
707 
708 }
mfem::IntegrationRule
Class for integration rule.
Definition: intrules.hpp:63
mfem::Geometry::NumGeom
static const int NumGeom
Definition: geom.hpp:34
mfem::Geometry::JacToPerfJac
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
Definition: geom.cpp:250
mfem::Mult
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:351
mfem::ElementTransformation::SetIntPoint
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:33
mfem::GeometryRefiner::Refine
RefinedGeometry * Refine(int Geom, int Times, int ETimes=1)
Definition: geom.cpp:289
mfem::DenseMatrix
Data type dense matrix.
Definition: densemat.hpp:22
mfem::Geometry::Volume
static const double Volume[NumGeom]
Definition: geom.hpp:37
mfem::Poly_1D::ClosedPoints
const double * ClosedPoints(const int p)
Definition: fe.cpp:6198
mfem::Geometry::NumBdrArray
static const int NumBdrArray[]
Definition: geom.hpp:35
mfem::RefinedGeometry::RefEdges
Array< int > RefEdges
Definition: geom.hpp:70
mfem::Geometry::GetVertices
const IntegrationRule * GetVertices(int GeomType)
Definition: geom.cpp:172
mfem::Geometries
Geometry Geometries
Definition: geom.cpp:267
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:205
mfem::GeometryRefiner::RefineInterior
const IntegrationRule * RefineInterior(int Geom, int Times)
Definition: geom.cpp:634
mfem::GlobGeometryRefiner
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:706
mfem::Linear3DFiniteElement
Class for linear FE on tetrahedron.
Definition: fe.hpp:643
mfem::RefinedGeometry::RefPts
IntegrationRule RefPts
Definition: geom.hpp:69
mfem::Linear2DFiniteElement
Class for linear FE on triangle.
Definition: fe.hpp:362
mfem::Geometry::Name
static const char * Name[NumGeom]
Definition: geom.hpp:36
mfem::IsoparametricTransformation::Jacobian
virtual const DenseMatrix & Jacobian()
Definition: eltrans.cpp:50
mfem::CalcInverse
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2588
mfem::Geometry::GetPerfPointMat
void GetPerfPointMat(int GeomType, DenseMatrix &pm)
Definition: geom.cpp:189
mfem::mfem_error
void mfem_error(const char *msg)
Definition: error.cpp:23
mfem::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
mfem::IsoparametricTransformation::SetFE
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:67
mfem::DenseMatrix::SetSize
void SetSize(int s)
If the matrix is not a square matrix of size s then recreate it.
Definition: densemat.cpp:73
mfem::RefinedGeometry::RefGeoms
Array< int > RefGeoms
Definition: geom.hpp:70
mfem::poly1d
Poly_1D poly1d
Definition: fe.cpp:6269
Generated on Fri Feb 18 2022 18:57:40 for MFEM by   doxygen 1.8.5