MFEM  v3.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mesh.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 // Implementation of data type mesh
13 
14 #include "mesh_headers.hpp"
15 #include "../fem/fem.hpp"
16 #include "../general/sort_pairs.hpp"
17 
18 #include <iostream>
19 #include <sstream>
20 #include <fstream>
21 #include <limits>
22 #include <cmath>
23 #include <cstring>
24 #include <ctime>
25 
26 namespace mfem
27 {
28 
29 using namespace std;
30 
32 {
33  int geom = GetElementBaseGeometry(i);
34  ElementTransformation *eltransf = GetElementTransformation(i);
35  eltransf->SetIntPoint(&Geometries.GetCenter(geom));
36  Geometries.JacToPerfJac(geom, eltransf->Jacobian(), J);
37 }
38 
39 double Mesh::GetElementSize(int i, int type)
40 {
41  DenseMatrix J(Dim);
42  GetElementJacobian(i, J);
43  if (type == 0)
44  return pow(fabs(J.Det()), 1./Dim);
45  else if (type == 1)
46  return J.CalcSingularvalue(Dim-1); // h_min
47  else
48  return J.CalcSingularvalue(0); // h_max
49 }
50 
51 double Mesh::GetElementSize(int i, const Vector &dir)
52 {
53  DenseMatrix J(Dim);
54  Vector d_hat(Dim);
55  GetElementJacobian(i, J);
56  J.MultTranspose(dir, d_hat);
57  return sqrt((d_hat * d_hat) / (dir * dir));
58 }
59 
61 {
62  ElementTransformation *et = GetElementTransformation(i);
63  const IntegrationRule &ir = IntRules.Get(GetElementBaseGeometry(i),
64  et->OrderJ());
65  double volume = 0.0;
66  for (int j = 0; j < ir.GetNPoints(); j++)
67  {
68  const IntegrationPoint &ip = ir.IntPoint(j);
69  et->SetIntPoint(&ip);
70  volume += ip.weight * et->Weight();
71  }
72 
73  return volume;
74 }
75 
77 {
78  int i, dim;
79  DenseMatrix J;
80  double h_min, h_max, kappa_min, kappa_max, h, kappa;
81 
82  cout << "Mesh Characteristics:" << flush;
83 
84  dim = Dimension();
85  J.SetSize(dim);
86 
87  if (Vh) Vh->SetSize(NumOfElements);
88  if (Vk) Vk->SetSize(NumOfElements);
89 
90  h_min = kappa_min = numeric_limits<double>::infinity();
91  h_max = kappa_max = -h_min;
92  for (i = 0; i < NumOfElements; i++)
93  {
94  GetElementJacobian(i, J);
95  h = pow(fabs(J.Det()), 1.0/double(dim));
96  kappa = J.CalcSingularvalue(0) / J.CalcSingularvalue(dim-1);
97  if (Vh) (*Vh)(i) = h;
98  if (Vk) (*Vk)(i) = kappa;
99 
100  if (h < h_min) h_min = h;
101  if (h > h_max) h_max = h;
102  if (kappa < kappa_min) kappa_min = kappa;
103  if (kappa > kappa_max) kappa_max = kappa;
104  }
105 
106  if (dim == 1)
107  cout << endl
108  << "Number of vertices : " << GetNV() << endl
109  << "Number of elements : " << GetNE() << endl
110  << "Number of bdr elem : " << GetNBE() << endl
111  << "h_min : " << h_min << endl
112  << "h_max : " << h_max << endl
113  << endl;
114  else if (dim == 2)
115  cout << endl
116  << "Number of vertices : " << GetNV() << endl
117  << "Number of edges : " << GetNEdges() << endl
118  << "Number of elements : " << GetNE() << endl
119  << "Number of bdr elem : " << GetNBE() << endl
120  << "Euler Number : " << EulerNumber2D() << endl
121  << "h_min : " << h_min << endl
122  << "h_max : " << h_max << endl
123  << "kappa_min : " << kappa_min << endl
124  << "kappa_max : " << kappa_max << endl
125  << endl;
126  else
127  cout << endl
128  << "Number of vertices : " << GetNV() << endl
129  << "Number of edges : " << GetNEdges() << endl
130  << "Number of faces : " << GetNFaces() << endl
131  << "Number of elements : " << GetNE() << endl
132  << "Number of bdr elem : " << GetNBE() << endl
133  << "Euler Number : " << EulerNumber() << endl
134  << "h_min : " << h_min << endl
135  << "h_max : " << h_max << endl
136  << "kappa_min : " << kappa_min << endl
137  << "kappa_max : " << kappa_max << endl
138  << endl;
139 }
140 
142 {
143  switch (ElemType)
144  {
145  case Element::POINT : return &PointFE;
146  case Element::SEGMENT : return &SegmentFE;
147  case Element::TRIANGLE : return &TriangleFE;
149  case Element::TETRAHEDRON : return &TetrahedronFE;
150  case Element::HEXAHEDRON : return &HexahedronFE;
151  }
152  mfem_error("Mesh::GetTransformationFEforElement - unknown ElementType");
153  return &TriangleFE;
154 }
155 
156 
158 {
159  ElTr->Attribute = GetAttribute(i);
160  ElTr->ElementNo = i;
161  if (Nodes == NULL)
162  {
163  GetPointMatrix(i, ElTr->GetPointMat());
164  ElTr->SetFE(GetTransformationFEforElementType(GetElementType(i)));
165  }
166  else
167  {
168  DenseMatrix &pm = ElTr->GetPointMat();
169  Array<int> vdofs;
170  Nodes->FESpace()->GetElementVDofs(i, vdofs);
171 
172  int n = vdofs.Size()/spaceDim;
173  pm.SetSize(spaceDim, n);
174  for (int k = 0; k < spaceDim; k++)
175  for (int j = 0; j < n; j++)
176  pm(k,j) = (*Nodes)(vdofs[n*k+j]);
177  ElTr->SetFE(Nodes->FESpace()->GetFE(i));
178  }
179 }
180 
181 void Mesh::GetElementTransformation(int i, const Vector &nodes,
183 {
184  ElTr->Attribute = GetAttribute(i);
185  ElTr->ElementNo = i;
186  DenseMatrix &pm = ElTr->GetPointMat();
187  if (Nodes == NULL)
188  {
189  int nv = elements[i]->GetNVertices();
190  const int *v = elements[i]->GetVertices();
191  int n = vertices.Size();
192  pm.SetSize(spaceDim, nv);
193  for (int k = 0; k < spaceDim; k++)
194  for (int j = 0; j < nv; j++)
195  pm(k, j) = nodes(k*n+v[j]);
196  ElTr->SetFE(GetTransformationFEforElementType(GetElementType(i)));
197  }
198  else
199  {
200  Array<int> vdofs;
201  Nodes->FESpace()->GetElementVDofs(i, vdofs);
202  int n = vdofs.Size()/spaceDim;
203  pm.SetSize(spaceDim, n);
204  for (int k = 0; k < spaceDim; k++)
205  for (int j = 0; j < n; j++)
206  pm(k,j) = nodes(vdofs[n*k+j]);
207  ElTr->SetFE(Nodes->FESpace()->GetFE(i));
208  }
209 }
210 
212 {
213  GetElementTransformation(i, &Transformation);
214 
215  return &Transformation;
216 }
217 
219 {
220  GetBdrElementTransformation(i, &FaceTransformation);
221  return &FaceTransformation;
222 }
223 
225 {
226  ElTr->Attribute = GetBdrAttribute(i);
227  ElTr->ElementNo = i; // boundary element number
228  if (Nodes == NULL)
229  {
230  GetBdrPointMatrix(i, ElTr->GetPointMat());
231  ElTr->SetFE(
232  GetTransformationFEforElementType(GetBdrElementType(i)));
233  }
234  else
235  {
236  DenseMatrix &pm = ElTr->GetPointMat();
237  Array<int> vdofs;
238  Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
239  int n = vdofs.Size()/spaceDim;
240  pm.SetSize(spaceDim, n);
241  for (int k = 0; k < spaceDim; k++)
242  for (int j = 0; j < n; j++)
243  pm(k,j) = (*Nodes)(vdofs[n*k+j]);
244  ElTr->SetFE(Nodes->FESpace()->GetBE(i));
245  }
246 }
247 
249 {
250  FTr->Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
251  FTr->ElementNo = FaceNo;
252  DenseMatrix &pm = FTr->GetPointMat();
253  if (Nodes == NULL)
254  {
255  const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
256  const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
257  pm.SetSize(spaceDim, nv);
258  for (int i = 0; i < spaceDim; i++)
259  for (int j = 0; j < nv; j++)
260  pm(i, j) = vertices[v[j]](i);
261  FTr->SetFE(GetTransformationFEforElementType(
262  (Dim == 1) ? Element::POINT : faces[FaceNo]->GetType()));
263  }
264  else
265  {
266  const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
267  if (face_el)
268  {
269  Array<int> vdofs;
270  Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
271  int n = vdofs.Size()/spaceDim;
272  pm.SetSize(spaceDim, n);
273  for (int i = 0; i < spaceDim; i++)
274  for (int j = 0; j < n; j++)
275  pm(i, j) = (*Nodes)(vdofs[n*i+j]);
276  FTr->SetFE(face_el);
277  }
278  else
279  {
280  int face_geom =
281  (Dim == 1) ? Geometry::POINT : faces[FaceNo]->GetGeometryType();
282  FaceInfo &face_info = faces_info[FaceNo];
283 
284  face_el = Nodes->FESpace()->GetTraceElement(face_info.Elem1No,
285  (Geometry::Type)face_geom);
286 
287  switch (face_geom)
288  {
289  case Geometry::POINT:
290  GetLocalPtToSegTransformation(FaceElemTr.Loc1.Transf,
291  face_info.Elem1Inf);
292  break;
293  case Geometry::SEGMENT:
294  if (GetElementType(face_info.Elem1No) == Element::TRIANGLE)
295  GetLocalSegToTriTransformation(FaceElemTr.Loc1.Transf,
296  face_info.Elem1Inf);
297  else // assume the element is a quad
298  GetLocalSegToQuadTransformation(FaceElemTr.Loc1.Transf,
299  face_info.Elem1Inf);
300  break;
301  case Geometry::TRIANGLE:
302  // --- assume the face is a triangle -- face of a tetrahedron
303  GetLocalTriToTetTransformation(FaceElemTr.Loc1.Transf,
304  face_info.Elem1Inf);
305  break;
306  case Geometry::SQUARE:
307  // --- assume the face is a quad -- face of a hexahedron
308  GetLocalQuadToHexTransformation(FaceElemTr.Loc1.Transf,
309  face_info.Elem1Inf);
310  break;
311  }
312 
313  IntegrationRule eir(face_el->GetDof());
314  FaceElemTr.Loc1.Transform(face_el->GetNodes(), eir);
315  // 'Transformation' is not used
316  Nodes->GetVectorValues(Transformation, eir, pm);
317 
318  FTr->SetFE(face_el);
319  }
320  }
321 }
322 
324 {
325  GetFaceTransformation(FaceNo, &FaceTransformation);
326  return &FaceTransformation;
327 }
328 
330 {
331  if (Dim == 2)
332  {
333  GetFaceTransformation(EdgeNo, EdTr);
334  return;
335  }
336  if (Dim == 1)
337  mfem_error("Mesh::GetEdgeTransformation not defined in 1D \n");
338 
339  EdTr->Attribute = 1;
340  EdTr->ElementNo = EdgeNo;
341  DenseMatrix &pm = EdTr->GetPointMat();
342  if (Nodes == NULL)
343  {
344  Array<int> v;
345  GetEdgeVertices(EdgeNo, v);
346  const int nv = 2;
347  pm.SetSize(spaceDim, nv);
348  for (int i = 0; i < spaceDim; i++)
349  for (int j = 0; j < nv; j++)
350  pm(i, j) = vertices[v[j]](i);
351  EdTr->SetFE(GetTransformationFEforElementType(Element::SEGMENT));
352 
353  }
354  else
355  {
356  Array<int> vdofs;
357  Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
358  int n = vdofs.Size()/spaceDim;
359  pm.SetSize(spaceDim, n);
360  for (int i = 0; i < spaceDim; i++)
361  for (int j = 0; j < n; j++)
362  pm(i, j) = (*Nodes)(vdofs[n*i+j]);
363  EdTr->SetFE(GetTransformationFEforElementType(Element::SEGMENT));
364  }
365 }
366 
368 {
369  GetEdgeTransformation(EdgeNo, &EdgeTransformation);
370  return &EdgeTransformation;
371 }
372 
373 
375  IsoparametricTransformation &Transf, int i)
376 {
377  const IntegrationRule *SegVert;
378  DenseMatrix &locpm = Transf.GetPointMat();
379 
380  Transf.SetFE(&PointFE);
382  locpm.SetSize(1, 1);
383  locpm(0, 0) = SegVert->IntPoint(i/64).x;
384  // (i/64) is the local face no. in the segment
385  // (i%64) is the orientation of the point (not used)
386 }
387 
389  IsoparametricTransformation &Transf, int i)
390 {
391  // tri_faces is the same as Triangle::edges
392  static const int tri_faces[3][2] = {{0, 1}, {1, 2}, {2, 0}};
393  static const int seg_inv_orient[2][2] = {{0, 1}, {1, 0}};
394  int j;
395  const int *tv, *so;
396  const IntegrationRule *TriVert;
397  DenseMatrix &locpm = Transf.GetPointMat();
398 
399  Transf.SetFE(&SegmentFE);
400  tv = tri_faces[i/64]; // (i/64) is the local face no. in the triangle
401  so = seg_inv_orient[i%64]; // (i%64) is the orientation of the segment
403  locpm.SetSize(2, 2);
404  for (j = 0; j < 2; j++)
405  {
406  locpm(0, so[j]) = TriVert->IntPoint(tv[j]).x;
407  locpm(1, so[j]) = TriVert->IntPoint(tv[j]).y;
408  }
409 }
410 
412  IsoparametricTransformation &Transf, int i)
413 {
414  // quad_faces is the same as Quadrilateral::edges
415  static const int quad_faces[4][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}};
416  static const int seg_inv_orient[2][2] = {{0, 1}, {1, 0}};
417  int j;
418  const int *qv, *so;
419  const IntegrationRule *QuadVert;
420  DenseMatrix &locpm = Transf.GetPointMat();
421 
422  Transf.SetFE(&SegmentFE);
423  qv = quad_faces[i/64]; // (i/64) is the local face no. in the quad
424  so = seg_inv_orient[i%64]; // (i%64) is the orientation of the segment
426  locpm.SetSize(2, 2);
427  for (j = 0; j < 2; j++)
428  {
429  locpm(0, so[j]) = QuadVert->IntPoint(qv[j]).x;
430  locpm(1, so[j]) = QuadVert->IntPoint(qv[j]).y;
431  }
432 }
433 
434 const int Mesh::tet_faces[4][3] =
435 {{1, 2, 3}, {0, 3, 2},
436  {0, 1, 3}, {0, 2, 1}};
437 
438 // same as Hexahedron::faces
439 const int Mesh::hex_faces[6][4] =
440 {{3, 2, 1, 0}, {0, 1, 5, 4},
441  {1, 2, 6, 5}, {2, 3, 7, 6},
442  {3, 0, 4, 7}, {4, 5, 6, 7}};
443 
444 const int Mesh::tri_orientations[6][3] =
445 {{0, 1, 2}, {1, 0, 2},
446  {2, 0, 1}, {2, 1, 0},
447  {1, 2, 0}, {0, 2, 1}};
448 
449 const int Mesh::quad_orientations[8][4] =
450 {{0, 1, 2, 3}, {0, 3, 2, 1},
451  {1, 2, 3, 0}, {1, 0, 3, 2},
452  {2, 3, 0, 1}, {2, 1, 0, 3},
453  {3, 0, 1, 2}, {3, 2, 1, 0}};
454 
456  IsoparametricTransformation &Transf, int i)
457 {
458  DenseMatrix &locpm = Transf.GetPointMat();
459 
460  Transf.SetFE(&TriangleFE);
461  // (i/64) is the local face no. in the tet
462  const int *tv = tet_faces[i/64];
463  // (i%64) is the orientation of the tetrahedron face
464  // w.r.t. the face element
465  const int *to = tri_orientations[i%64];
466  const IntegrationRule *TetVert =
468  locpm.SetSize(3, 3);
469  for (int j = 0; j < 3; j++)
470  {
471  const IntegrationPoint &vert = TetVert->IntPoint(tv[to[j]]);
472  locpm(0, j) = vert.x;
473  locpm(1, j) = vert.y;
474  locpm(2, j) = vert.z;
475  }
476 }
477 
479  IsoparametricTransformation &Transf, int i)
480 {
481  DenseMatrix &locpm = Transf.GetPointMat();
482 
483  Transf.SetFE(&QuadrilateralFE);
484  // (i/64) is the local face no. in the hex
485  const int *hv = hex_faces[i/64];
486  // (i%64) is the orientation of the quad
487  const int *qo = quad_orientations[i%64];
489  locpm.SetSize(3, 4);
490  for (int j = 0; j < 4; j++)
491  {
492  const IntegrationPoint &vert = HexVert->IntPoint(hv[qo[j]]);
493  locpm(0, j) = vert.x;
494  locpm(1, j) = vert.y;
495  locpm(2, j) = vert.z;
496  }
497 }
498 
500  int mask)
501 {
502  FaceInfo &face_info = faces_info[FaceNo];
503 
504  // setup the transformation for the first element
505  FaceElemTr.Elem1No = face_info.Elem1No;
506  if (mask & 1)
507  {
508  GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
509  FaceElemTr.Elem1 = &Transformation;
510  }
511  else
512  FaceElemTr.Elem1 = NULL;
513 
514  // setup the transformation for the second element
515  // return NULL in the Elem2 field if there's no second element, i.e.
516  // the face is on the "boundary"
517  if ( (FaceElemTr.Elem2No = face_info.Elem2No) >= 0 && (mask & 2))
518  {
519 #ifdef MFEM_DEBUG
520  if (NURBSext && (mask & 1))
521  mfem_error("Mesh::GetFaceElementTransformations :"
522  " NURBS mesh is not supported!");
523 #endif
524  GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
525  FaceElemTr.Elem2 = &Transformation2;
526  }
527  else
528  FaceElemTr.Elem2 = NULL;
529 
530  if (Dim == 1)
531  FaceElemTr.FaceGeom = Geometry::POINT;
532  else
533  FaceElemTr.FaceGeom = faces[FaceNo]->GetGeometryType();
534 
535  // setup the face transformation
536  if (mask & 16)
537  FaceElemTr.Face = GetFaceTransformation(FaceNo);
538  else
539  FaceElemTr.Face = NULL;
540 
541  // setup Loc1 & Loc2
542  int face_type = (Dim == 1) ? Element::POINT : faces[FaceNo]->GetType();
543  switch (face_type)
544  {
545  case Element::POINT:
546  if (mask & 4)
547  GetLocalPtToSegTransformation(FaceElemTr.Loc1.Transf,
548  face_info.Elem1Inf);
549 
550  if (FaceElemTr.Elem2No >= 0 && (mask & 8))
551  GetLocalPtToSegTransformation(FaceElemTr.Loc2.Transf,
552  face_info.Elem2Inf);
553  break;
554  case Element::SEGMENT:
555  if (mask & 4)
556  {
557  if (GetElementType(face_info.Elem1No) == Element::TRIANGLE)
558  GetLocalSegToTriTransformation(FaceElemTr.Loc1.Transf,
559  face_info.Elem1Inf);
560  else // assume the element is a quad
561  GetLocalSegToQuadTransformation(FaceElemTr.Loc1.Transf,
562  face_info.Elem1Inf);
563  }
564 
565  if (FaceElemTr.Elem2No >= 0 && (mask & 8))
566  {
567  if (GetElementType(face_info.Elem2No) == Element::TRIANGLE)
568  GetLocalSegToTriTransformation(FaceElemTr.Loc2.Transf,
569  face_info.Elem2Inf);
570  else // assume the element is a quad
571  GetLocalSegToQuadTransformation(FaceElemTr.Loc2.Transf,
572  face_info.Elem2Inf);
573  }
574  break;
575  case Element::TRIANGLE:
576  // --------- assumes the face is a triangle -- face of a tetrahedron
577  if (mask & 4)
578  GetLocalTriToTetTransformation(FaceElemTr.Loc1.Transf,
579  face_info.Elem1Inf);
580  if (FaceElemTr.Elem2No >= 0 && (mask & 8))
581  GetLocalTriToTetTransformation(FaceElemTr.Loc2.Transf,
582  face_info.Elem2Inf);
583  break;
585  // --------- assumes the face is a quad -- face of a hexahedron
586  if (mask & 4)
587  GetLocalQuadToHexTransformation(FaceElemTr.Loc1.Transf,
588  face_info.Elem1Inf);
589  if (FaceElemTr.Elem2No >= 0 && (mask & 8))
590  GetLocalQuadToHexTransformation(FaceElemTr.Loc2.Transf,
591  face_info.Elem2Inf);
592  break;
593  }
594 
595  return &FaceElemTr;
596 }
597 
599 {
601  int fn;
602  if (Dim == 3)
603  fn = be_to_face[BdrElemNo];
604  else if (Dim == 2)
605  fn = be_to_edge[BdrElemNo];
606  else
607  fn = boundary[BdrElemNo]->GetVertices()[0];
608  if (FaceIsTrueInterior(fn))
609  return NULL;
610  tr = GetFaceElementTransformations(fn);
611  tr->Face->Attribute = boundary[BdrElemNo]->GetAttribute();
612  return tr;
613 }
614 
615 void Mesh::GetFaceElements(int Face, int *Elem1, int *Elem2)
616 {
617  *Elem1 = faces_info[Face].Elem1No;
618  *Elem2 = faces_info[Face].Elem2No;
619 }
620 
621 void Mesh::GetFaceInfos(int Face, int *Inf1, int *Inf2)
622 {
623  *Inf1 = faces_info[Face].Elem1Inf;
624  *Inf2 = faces_info[Face].Elem2Inf;
625 }
626 
628 {
629  NumOfVertices = NumOfElements = NumOfBdrElements = NumOfEdges = -1;
630  WantTwoLevelState = 0;
631  State = Mesh::NORMAL;
632  Nodes = NULL;
633  own_nodes = 1;
634  NURBSext = NULL;
635  ncmesh = NULL;
636  nc_coarse_level = NULL;
637 }
638 
640 {
641  el_to_edge = el_to_face = el_to_el =
642  bel_to_edge = face_edge = edge_vertex = NULL;
643 }
644 
646 {
647  if (el_to_edge != NULL)
648  delete el_to_edge;
649 
650  if (el_to_face != NULL)
651  delete el_to_face;
652 
653  if (el_to_el != NULL)
654  delete el_to_el;
655 
656  if (Dim == 3 && bel_to_edge != NULL)
657  delete bel_to_edge;
658 
659  if (face_edge != NULL)
660  delete face_edge;
661 
662  if (edge_vertex != NULL)
663  delete edge_vertex;
664 
665  InitTables();
666 
667  delete nc_coarse_level;
668  nc_coarse_level = NULL;
669 }
670 
672 {
673  delete el_to_el;
674  delete face_edge;
675  delete edge_vertex;
676 
677  el_to_el = face_edge = edge_vertex = NULL;
678 
679  delete nc_coarse_level;
680  nc_coarse_level = NULL;
681 }
682 
684 {
685  int i, j, nattr;
686  Array<int> attribs;
687 
688  attribs.SetSize(GetNBE());
689  for (i = 0; i < attribs.Size(); i++)
690  attribs[i] = GetBdrAttribute(i);
691  attribs.Sort();
692 
693  if (attribs.Size() > 0)
694  nattr = 1;
695  else
696  nattr = 0;
697  for (i = 1; i < attribs.Size(); i++)
698  if (attribs[i] != attribs[i-1])
699  nattr++;
700 
701  bdr_attributes.SetSize(nattr);
702  if (nattr > 0)
703  {
704  bdr_attributes[0] = attribs[0];
705  for (i = j = 1; i < attribs.Size(); i++)
706  if (attribs[i] != attribs[i-1])
707  bdr_attributes[j++] = attribs[i];
708  if (attribs[0] <= 0)
709  cout << "Mesh::SetAttributes(): "
710  "Non-positive attributes on the boundary!"
711  << endl;
712  }
713 
714 
715  attribs.SetSize(GetNE());
716  for (i = 0; i < attribs.Size(); i++)
717  attribs[i] = GetAttribute(i);
718  attribs.Sort();
719 
720  if (attribs.Size() > 0)
721  nattr = 1;
722  else
723  nattr = 0;
724  for (i = 1; i < attribs.Size(); i++)
725  if (attribs[i] != attribs[i-1])
726  nattr++;
727 
728  attributes.SetSize(nattr);
729  if (nattr > 0)
730  {
731  attributes[0] = attribs[0];
732  for (i = j = 1; i < attribs.Size(); i++)
733  if (attribs[i] != attribs[i-1])
734  attributes[j++] = attribs[i];
735  if (attribs[0] <= 0)
736  cout << "Mesh::SetAttributes(): "
737  "Non-positive attributes in the domain!"
738  << endl;
739  }
740 }
741 
742 void Mesh::InitMesh(int _Dim, int _spaceDim, int NVert, int NElem, int NBdrElem)
743 {
744  Dim = _Dim;
745  spaceDim = _spaceDim;
746 
747  Init();
748  InitTables();
749 
750  NumOfVertices = 0;
751  vertices.SetSize(NVert); // just allocate space for vertices
752 
753  NumOfElements = 0;
754  elements.SetSize(NElem); // just allocate space for Element *
755 
756  NumOfBdrElements = 0;
757  boundary.SetSize(NBdrElem); // just allocate space for Element *
758 }
759 
760 void Mesh::AddVertex(const double *x)
761 {
762  double *y = vertices[NumOfVertices]();
763 
764  for (int i = 0; i < spaceDim; i++)
765  y[i] = x[i];
766  NumOfVertices++;
767 }
768 
769 void Mesh::AddTri(const int *vi, int attr)
770 {
771  elements[NumOfElements++] = new Triangle(vi, attr);
772 }
773 
774 void Mesh::AddTriangle(const int *vi, int attr)
775 {
776  elements[NumOfElements++] = new Triangle(vi, attr);
777 }
778 
779 void Mesh::AddQuad(const int *vi, int attr)
780 {
781  elements[NumOfElements++] = new Quadrilateral(vi, attr);
782 }
783 
784 void Mesh::AddTet(const int *vi, int attr)
785 {
786 #ifdef MFEM_USE_MEMALLOC
787  Tetrahedron *tet;
788  tet = TetMemory.Alloc();
789  tet->SetVertices(vi);
790  tet->SetAttribute(attr);
791  elements[NumOfElements++] = tet;
792 #else
793  elements[NumOfElements++] = new Tetrahedron(vi, attr);
794 #endif
795 }
796 
797 void Mesh::AddHex(const int *vi, int attr)
798 {
799  elements[NumOfElements++] = new Hexahedron(vi, attr);
800 }
801 
802 void Mesh::AddHexAsTets(const int *vi, int attr)
803 {
804  static const int hex_to_tet[6][4] =
805  { { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
806  { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 } };
807  int ti[4];
808 
809  for (int i = 0; i < 6; i++)
810  {
811  for (int j = 0; j < 4; j++)
812  ti[j] = vi[hex_to_tet[i][j]];
813  AddTet(ti, attr);
814  }
815 }
816 
817 void Mesh::AddBdrSegment(const int *vi, int attr)
818 {
819  boundary[NumOfBdrElements++] = new Segment(vi, attr);
820 }
821 
822 void Mesh::AddBdrTriangle(const int *vi, int attr)
823 {
824  boundary[NumOfBdrElements++] = new Triangle(vi, attr);
825 }
826 
827 void Mesh::AddBdrQuad(const int *vi, int attr)
828 {
829  boundary[NumOfBdrElements++] = new Quadrilateral(vi, attr);
830 }
831 
832 void Mesh::AddBdrQuadAsTriangles(const int *vi, int attr)
833 {
834  static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
835  int ti[3];
836 
837  for (int i = 0; i < 2; i++)
838  {
839  for (int j = 0; j < 3; j++)
840  ti[j] = vi[quad_to_tri[i][j]];
841  AddBdrTriangle(ti, attr);
842  }
843 }
844 
846 {
847  int i, j;
848  Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
849 
850  // GenerateFaces();
851 
852  for (i = 0; i < boundary.Size(); i++)
853  FreeElement(boundary[i]);
854 
855  if (Dim == 3)
856  {
857  delete bel_to_edge;
858  bel_to_edge = NULL;
859  }
860 
861  // count the 'NumOfBdrElements'
862  NumOfBdrElements = 0;
863  for (i = 0; i < faces_info.Size(); i++)
864  if (faces_info[i].Elem2No < 0)
865  NumOfBdrElements++;
866 
867  boundary.SetSize(NumOfBdrElements);
868  be2face.SetSize(NumOfBdrElements);
869  for (j = i = 0; i < faces_info.Size(); i++)
870  if (faces_info[i].Elem2No < 0)
871  {
872  boundary[j] = faces[i]->Duplicate(this);
873  be2face[j++] = i;
874  }
875  // In 3D, 'bel_to_edge' is destroyed but it's not updated.
876 }
877 
878 typedef struct {
879  int edge;
880  double length;
881 } edge_length;
882 
883 // Used by qsort to sort edges in increasing (according their length) order.
884 static int edge_compare(const void *ii, const void *jj)
885 {
886  edge_length *i = (edge_length *)ii, *j = (edge_length *)jj;
887  if (i->length > j->length) return (1);
888  if (i->length < j->length) return (-1);
889  return (0);
890 }
891 
892 void Mesh::FinalizeTriMesh(int generate_edges, int refine, bool fix_orientation)
893 {
894  CheckElementOrientation(fix_orientation);
895 
896  if (refine)
897  MarkTriMeshForRefinement();
898 
899  if (generate_edges)
900  {
901  el_to_edge = new Table;
902  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
903  GenerateFaces();
904  CheckBdrElementOrientation();
905  }
906  else
907  NumOfEdges = 0;
908 
909  NumOfFaces = 0;
910 
911  SetAttributes();
912 
913  meshgen = 1;
914 }
915 
916 void Mesh::FinalizeQuadMesh(int generate_edges, int refine,
917  bool fix_orientation)
918 {
919  if (fix_orientation)
920  CheckElementOrientation(fix_orientation);
921 
922  if (generate_edges)
923  {
924  el_to_edge = new Table;
925  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
926  GenerateFaces();
927  CheckBdrElementOrientation();
928  }
929  else
930  NumOfEdges = 0;
931 
932  NumOfFaces = 0;
933 
934  SetAttributes();
935 
936  meshgen = 2;
937 }
938 
940 {
941  if (meshgen & 1)
942  {
943  if (Dim == 2)
944  MarkTriMeshForRefinement();
945  else if (Dim == 3)
946  MarkTetMeshForRefinement();
947  }
948 }
949 
951 {
952  // Mark the longest triangle edge by rotating the indeces so that
953  // vertex 0 - vertex 1 to be the longest element's edge.
954  DenseMatrix pmat;
955  for (int i = 0; i < NumOfElements; i++)
956  if (elements[i]->GetType() == Element::TRIANGLE)
957  {
958  GetPointMatrix(i, pmat);
959  elements[i]->MarkEdge(pmat);
960  }
961 }
962 
964 {
965  NumOfEdges = v_to_v.NumberOfEntries();
966  edge_length *length = new edge_length[NumOfEdges];
967  for (int i = 0; i < NumOfVertices; i++)
968  {
969  for (DSTable::RowIterator it(v_to_v, i); !it; ++it)
970  {
971  int j = it.Index();
972  length[j].length = GetLength(i, it.Column());
973  length[j].edge = j;
974  }
975  }
976 
977  // sort in increasing order
978  qsort(length, NumOfEdges, sizeof(edge_length), edge_compare);
979 
980  order.SetSize(NumOfEdges);
981  for (int i = 0; i < NumOfEdges; i++)
982  order[length[i].edge] = i;
983 
984  delete [] length;
985 }
986 
988 {
989  // Mark the longest tetrahedral edge by rotating the indices so that
990  // vertex 0 - vertex 1 is the longest edge in the element.
991  DSTable v_to_v(NumOfVertices);
992  GetVertexToVertexTable(v_to_v);
993  Array<int> order;
994 
995  GetEdgeOrdering(v_to_v, order);
996 
997  for (int i = 0; i < NumOfElements; i++)
998  if (elements[i]->GetType() == Element::TETRAHEDRON)
999  elements[i]->MarkEdge(v_to_v, order);
1000 
1001  for (int i = 0; i < NumOfBdrElements; i++)
1002  if (boundary[i]->GetType() == Element::TRIANGLE)
1003  boundary[i]->MarkEdge(v_to_v, order);
1004 }
1005 
1006 void Mesh::PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert)
1007 {
1008  if (*old_v_to_v && *old_elem_vert)
1009  return;
1010 
1011  FiniteElementSpace *fes = Nodes->FESpace();
1012  const FiniteElementCollection *fec = fes->FEColl();
1013 
1014  if (*old_v_to_v == NULL)
1015  {
1016  int num_edge_dofs = fec->DofForGeometry(Geometry::SEGMENT);
1017  if (num_edge_dofs > 0)
1018  {
1019  *old_v_to_v = new DSTable(NumOfVertices);
1020  GetVertexToVertexTable(*(*old_v_to_v));
1021  }
1022  }
1023  if (*old_elem_vert == NULL)
1024  {
1025  // assuming all elements have the same geometry
1026  int num_elem_dofs = fec->DofForGeometry(GetElementBaseGeometry(0));
1027  if (num_elem_dofs > 1)
1028  {
1029  *old_elem_vert = new Table;
1030  (*old_elem_vert)->MakeI(GetNE());
1031  for (int i = 0; i < GetNE(); i++)
1032  {
1033  (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
1034  }
1035  (*old_elem_vert)->MakeJ();
1036  for (int i = 0; i < GetNE(); i++)
1037  {
1038  (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
1039  elements[i]->GetNVertices());
1040  }
1041  (*old_elem_vert)->ShiftUpI();
1042  }
1043  }
1044 }
1045 
1046 void Mesh::DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
1047 {
1048  FiniteElementSpace *fes = Nodes->FESpace();
1049  const FiniteElementCollection *fec = fes->FEColl();
1050  int num_edge_dofs = fec->DofForGeometry(Geometry::SEGMENT);
1051  // assuming all faces have the same geometry
1052  int num_face_dofs =
1053  (Dim < 3) ? 0 : fec->DofForGeometry(GetFaceBaseGeometry(0));
1054  // assuming all elements have the same geometry
1055  int num_elem_dofs = fec->DofForGeometry(GetElementBaseGeometry(0));
1056 
1057  // reorder the Nodes
1058  Vector onodes = *Nodes;
1059 
1060  Array<int> old_dofs, new_dofs;
1061  int offset;
1062 #ifdef MFEM_DEBUG
1063  int redges = 0;
1064 #endif
1065 
1066  // vertex dofs do not need to be moved
1067  offset = NumOfVertices * fec->DofForGeometry(Geometry::POINT);
1068 
1069  // edge dofs:
1070  // edge enumeration may be different but edge orientation is the same
1071  if (num_edge_dofs > 0)
1072  {
1073  DSTable new_v_to_v(NumOfVertices);
1074  GetVertexToVertexTable(new_v_to_v);
1075 
1076  for (int i = 0; i < NumOfVertices; i++)
1077  {
1078  for (DSTable::RowIterator it(new_v_to_v, i); !it; ++it)
1079  {
1080  int old_i = (*old_v_to_v)(i, it.Column());
1081  int new_i = it.Index();
1082 #ifdef MFEM_DEBUG
1083  if (old_i != new_i)
1084  redges++;
1085 #endif
1086  old_dofs.SetSize(num_edge_dofs);
1087  new_dofs.SetSize(num_edge_dofs);
1088  for (int j = 0; j < num_edge_dofs; j++)
1089  {
1090  old_dofs[j] = offset + old_i * num_edge_dofs + j;
1091  new_dofs[j] = offset + new_i * num_edge_dofs + j;
1092  }
1093  fes->DofsToVDofs(old_dofs);
1094  fes->DofsToVDofs(new_dofs);
1095  for (int j = 0; j < old_dofs.Size(); j++)
1096  (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
1097  }
1098  }
1099  offset += NumOfEdges * num_edge_dofs;
1100  }
1101 #ifdef MFEM_DEBUG
1102  cout << "Mesh::DoNodeReorder : redges = " << redges << endl;
1103 #endif
1104 
1105  // face dofs:
1106  // both enumeration and orientation of the faces may be different
1107  if (num_face_dofs > 0)
1108  {
1109  // generate the old face-vertex table using the unmodified 'faces'
1110  Table old_face_vertex;
1111  old_face_vertex.MakeI(NumOfFaces);
1112  for (int i = 0; i < NumOfFaces; i++)
1113  old_face_vertex.AddColumnsInRow(i, faces[i]->GetNVertices());
1114  old_face_vertex.MakeJ();
1115  for (int i = 0; i < NumOfFaces; i++)
1116  old_face_vertex.AddConnections(i, faces[i]->GetVertices(),
1117  faces[i]->GetNVertices());
1118  old_face_vertex.ShiftUpI();
1119 
1120  // update 'el_to_face', 'be_to_face', 'faces', and 'faces_info'
1121  STable3D *faces_tbl = GetElementToFaceTable(1);
1122  GenerateFaces();
1123 
1124  // loop over the old face numbers
1125  for (int i = 0; i < NumOfFaces; i++)
1126  {
1127  int *old_v = old_face_vertex.GetRow(i), *new_v;
1128  int new_i, new_or, *dof_ord;
1129  switch (old_face_vertex.RowSize(i))
1130  {
1131  case 3:
1132  new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
1133  new_v = faces[new_i]->GetVertices();
1134  new_or = GetTriOrientation(old_v, new_v);
1135  dof_ord = fec->DofOrderForOrientation(Geometry::TRIANGLE, new_or);
1136  break;
1137  case 4:
1138  default:
1139  new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
1140  new_v = faces[new_i]->GetVertices();
1141  new_or = GetQuadOrientation(old_v, new_v);
1142  dof_ord = fec->DofOrderForOrientation(Geometry::SQUARE, new_or);
1143  break;
1144  }
1145 
1146  old_dofs.SetSize(num_face_dofs);
1147  new_dofs.SetSize(num_face_dofs);
1148  for (int j = 0; j < num_face_dofs; j++)
1149  {
1150  old_dofs[j] = offset + i * num_face_dofs + j;
1151  new_dofs[j] = offset + new_i * num_face_dofs + dof_ord[j];
1152  // we assumed the dofs are non-directional
1153  // i.e. dof_ord[j] is >= 0
1154  }
1155  fes->DofsToVDofs(old_dofs);
1156  fes->DofsToVDofs(new_dofs);
1157  for (int j = 0; j < old_dofs.Size(); j++)
1158  (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
1159  }
1160 
1161  offset += NumOfFaces * num_face_dofs;
1162  delete faces_tbl;
1163  }
1164 
1165  // element dofs:
1166  // element orientation may be different
1167  if (num_elem_dofs > 1)
1168  {
1169  // matters when the 'fec' is
1170  // (this code is executed only for triangles/tets)
1171  // - Pk on triangles, k >= 4
1172  // - Qk on quads, k >= 3
1173  // - Pk on tets, k >= 5
1174  // - Qk on hexes, k >= 3
1175  // - DG spaces
1176  // - ...
1177 
1178  // loop over all elements
1179  for (int i = 0; i < GetNE(); i++)
1180  {
1181  int *old_v = old_elem_vert->GetRow(i);
1182  int *new_v = elements[i]->GetVertices();
1183  int new_or, *dof_ord;
1184  int geom = elements[i]->GetGeometryType();
1185  switch (geom)
1186  {
1187  case Geometry::SEGMENT:
1188  new_or = (old_v[0] == new_v[0]) ? +1 : -1;
1189  break;
1190  case Geometry::TRIANGLE:
1191  new_or = GetTriOrientation(old_v, new_v);
1192  break;
1193  case Geometry::SQUARE:
1194  new_or = GetQuadOrientation(old_v, new_v);
1195  break;
1196  default:
1197  new_or = 0;
1198  cerr << "Mesh::DoNodeReorder : " << Geometry::Name[geom]
1199  << " elements (" << fec->Name()
1200  << " FE collection) are not supported yet!" << endl;
1201  mfem_error();
1202  break;
1203  }
1204  dof_ord = fec->DofOrderForOrientation(geom, new_or);
1205  if (dof_ord == NULL)
1206  {
1207  cerr << "Mesh::DoNodeReorder : FE collection '" << fec->Name()
1208  << "' does not define reordering for " << Geometry::Name[geom]
1209  << " elements!" << endl;
1210  mfem_error();
1211  }
1212  old_dofs.SetSize(num_elem_dofs);
1213  new_dofs.SetSize(num_elem_dofs);
1214  for (int j = 0; j < num_elem_dofs; j++)
1215  {
1216  // we assume the dofs are non-directional, i.e. dof_ord[j] is >= 0
1217  old_dofs[j] = offset + dof_ord[j];
1218  new_dofs[j] = offset + j;
1219  }
1220  fes->DofsToVDofs(old_dofs);
1221  fes->DofsToVDofs(new_dofs);
1222  for (int j = 0; j < old_dofs.Size(); j++)
1223  (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
1224 
1225  offset += num_elem_dofs;
1226  }
1227  }
1228 
1229  // Update Tables, faces, etc
1230  if (Dim > 2)
1231  {
1232  if (num_face_dofs == 0)
1233  {
1234  // needed for FE spaces that have face dofs, even if
1235  // the 'Nodes' do not have face dofs.
1236  GetElementToFaceTable();
1237  GenerateFaces();
1238  }
1239  CheckBdrElementOrientation();
1240  }
1241  if (el_to_edge)
1242  {
1243  // update 'el_to_edge', 'be_to_edge' (2D), 'bel_to_edge' (3D)
1244  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1245  if (Dim == 2)
1246  {
1247  // update 'faces' and 'faces_info'
1248  GenerateFaces();
1249  CheckBdrElementOrientation();
1250  }
1251  }
1252 }
1253 
1254 void Mesh::FinalizeTetMesh(int generate_edges, int refine, bool fix_orientation)
1255 {
1256  CheckElementOrientation(fix_orientation);
1257 
1258  if (NumOfBdrElements == 0)
1259  {
1260  GetElementToFaceTable();
1261  GenerateFaces();
1262  GenerateBoundaryElements();
1263  }
1264 
1265  if (refine)
1266  {
1267  MarkTetMeshForRefinement();
1268  }
1269 
1270  GetElementToFaceTable();
1271  GenerateFaces();
1272 
1273  CheckBdrElementOrientation();
1274 
1275  if (generate_edges == 1)
1276  {
1277  el_to_edge = new Table;
1278  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1279  }
1280  else
1281  {
1282  el_to_edge = NULL; // Not really necessary -- InitTables was called
1283  bel_to_edge = NULL;
1284  NumOfEdges = 0;
1285  }
1286 
1287  SetAttributes();
1288 
1289  meshgen = 1;
1290 }
1291 
1292 void Mesh::FinalizeHexMesh(int generate_edges, int refine, bool fix_orientation)
1293 {
1294  CheckElementOrientation(fix_orientation);
1295 
1296  GetElementToFaceTable();
1297  GenerateFaces();
1298 
1299  if (NumOfBdrElements == 0)
1300  GenerateBoundaryElements();
1301 
1302  CheckBdrElementOrientation();
1303 
1304  if (generate_edges)
1305  {
1306  el_to_edge = new Table;
1307  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1308  }
1309  else
1310  NumOfEdges = 0;
1311 
1312  SetAttributes();
1313 
1314  meshgen = 2;
1315 }
1316 
1317 void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
1318  int generate_edges, double sx, double sy, double sz)
1319 {
1320  int x, y, z;
1321 
1322  int NVert, NElem, NBdrElem;
1323 
1324  NVert = (nx+1) * (ny+1) * (nz+1);
1325  NElem = nx * ny * nz;
1326  NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
1327  if (type == Element::TETRAHEDRON)
1328  {
1329  NElem *= 6;
1330  NBdrElem *= 2;
1331  }
1332 
1333  InitMesh(3, 3, NVert, NElem, NBdrElem);
1334 
1335  double coord[3];
1336  int ind[8];
1337 
1338  // Sets vertices and the corresponding coordinates
1339  for (z = 0; z <= nz; z++)
1340  {
1341  coord[2] = ((double) z / nz) * sz;
1342  for (y = 0; y <= ny; y++)
1343  {
1344  coord[1] = ((double) y / ny) * sy;
1345  for (x = 0; x <= nx; x++)
1346  {
1347  coord[0] = ((double) x / nx) * sx;
1348  AddVertex(coord);
1349  }
1350  }
1351  }
1352 
1353 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
1354 
1355  // Sets elements and the corresponding indices of vertices
1356  for (z = 0; z < nz; z++)
1357  {
1358  for (y = 0; y < ny; y++)
1359  {
1360  for (x = 0; x < nx; x++)
1361  {
1362  ind[0] = VTX(x , y , z );
1363  ind[1] = VTX(x+1, y , z );
1364  ind[2] = VTX(x+1, y+1, z );
1365  ind[3] = VTX(x , y+1, z );
1366  ind[4] = VTX(x , y , z+1);
1367  ind[5] = VTX(x+1, y , z+1);
1368  ind[6] = VTX(x+1, y+1, z+1);
1369  ind[7] = VTX(x , y+1, z+1);
1370  if (type == Element::TETRAHEDRON)
1371  AddHexAsTets(ind, 1);
1372  else
1373  AddHex(ind, 1);
1374  }
1375  }
1376  }
1377 
1378  // Sets boundary elements and the corresponding indices of vertices
1379  // bottom, bdr. attribute 1
1380  for (y = 0; y < ny; y++)
1381  for (x = 0; x < nx; x++)
1382  {
1383  ind[0] = VTX(x , y , 0);
1384  ind[1] = VTX(x , y+1, 0);
1385  ind[2] = VTX(x+1, y+1, 0);
1386  ind[3] = VTX(x+1, y , 0);
1387  if (type == Element::TETRAHEDRON)
1388  AddBdrQuadAsTriangles(ind, 1);
1389  else
1390  AddBdrQuad(ind, 1);
1391  }
1392  // top, bdr. attribute 6
1393  for (y = 0; y < ny; y++)
1394  for (x = 0; x < nx; x++)
1395  {
1396  ind[0] = VTX(x , y , nz);
1397  ind[1] = VTX(x+1, y , nz);
1398  ind[2] = VTX(x+1, y+1, nz);
1399  ind[3] = VTX(x , y+1, nz);
1400  if (type == Element::TETRAHEDRON)
1401  AddBdrQuadAsTriangles(ind, 6);
1402  else
1403  AddBdrQuad(ind, 6);
1404  }
1405  // left, bdr. attribute 5
1406  for (z = 0; z < nz; z++)
1407  for (y = 0; y < ny; y++)
1408  {
1409  ind[0] = VTX(0 , y , z );
1410  ind[1] = VTX(0 , y , z+1);
1411  ind[2] = VTX(0 , y+1, z+1);
1412  ind[3] = VTX(0 , y+1, z );
1413  if (type == Element::TETRAHEDRON)
1414  AddBdrQuadAsTriangles(ind, 5);
1415  else
1416  AddBdrQuad(ind, 5);
1417  }
1418  // right, bdr. attribute 3
1419  for (z = 0; z < nz; z++)
1420  for (y = 0; y < ny; y++)
1421  {
1422  ind[0] = VTX(nx, y , z );
1423  ind[1] = VTX(nx, y+1, z );
1424  ind[2] = VTX(nx, y+1, z+1);
1425  ind[3] = VTX(nx, y , z+1);
1426  if (type == Element::TETRAHEDRON)
1427  AddBdrQuadAsTriangles(ind, 3);
1428  else
1429  AddBdrQuad(ind, 3);
1430  }
1431  // front, bdr. attribute 2
1432  for (x = 0; x < nx; x++)
1433  for (z = 0; z < nz; z++)
1434  {
1435  ind[0] = VTX(x , 0, z );
1436  ind[1] = VTX(x+1, 0, z );
1437  ind[2] = VTX(x+1, 0, z+1);
1438  ind[3] = VTX(x , 0, z+1);
1439  if (type == Element::TETRAHEDRON)
1440  AddBdrQuadAsTriangles(ind, 2);
1441  else
1442  AddBdrQuad(ind, 2);
1443  }
1444  // back, bdr. attribute 4
1445  for (x = 0; x < nx; x++)
1446  for (z = 0; z < nz; z++)
1447  {
1448  ind[0] = VTX(x , ny, z );
1449  ind[1] = VTX(x , ny, z+1);
1450  ind[2] = VTX(x+1, ny, z+1);
1451  ind[3] = VTX(x+1, ny, z );
1452  if (type == Element::TETRAHEDRON)
1453  AddBdrQuadAsTriangles(ind, 4);
1454  else
1455  AddBdrQuad(ind, 4);
1456  }
1457 
1458 #if 0
1459  ofstream test_stream("debug.mesh");
1460  Print(test_stream);
1461  test_stream.close();
1462 #endif
1463 
1464  int refine = 1;
1465  bool fix_orientation = true;
1466 
1467  if (type == Element::TETRAHEDRON)
1468  FinalizeTetMesh(generate_edges, refine, fix_orientation);
1469  else
1470  FinalizeHexMesh(generate_edges, refine, fix_orientation);
1471 }
1472 
1473 void Mesh::Make2D(int nx, int ny, Element::Type type, int generate_edges,
1474  double sx, double sy)
1475 {
1476  int i, j, k;
1477 
1478  Dim = spaceDim = 2;
1479 
1480  Init();
1481  InitTables();
1482 
1483  // Creates quadrilateral mesh
1484  if (type == Element::QUADRILATERAL)
1485  {
1486  meshgen = 2;
1487  NumOfVertices = (nx+1) * (ny+1);
1488  NumOfElements = nx * ny;
1489  NumOfBdrElements = 2 * nx + 2 * ny;
1490  vertices.SetSize(NumOfVertices);
1491  elements.SetSize(NumOfElements);
1492  boundary.SetSize(NumOfBdrElements);
1493  double cx, cy;
1494  int ind[4];
1495 
1496  // Sets vertices and the corresponding coordinates
1497  k = 0;
1498  for (j = 0; j < ny+1; j++)
1499  {
1500  cy = ((double) j / ny) * sy;
1501  for (i = 0; i < nx+1; i++)
1502  {
1503  cx = ((double) i / nx) * sx;
1504  vertices[k](0) = cx;
1505  vertices[k](1) = cy;
1506  k++;
1507  }
1508  }
1509 
1510  // Sets elements and the corresponding indices of vertices
1511  k = 0;
1512  for (j = 0; j < ny; j++)
1513  {
1514  for (i = 0; i < nx; i++)
1515  {
1516  ind[0] = i + j*(nx+1);
1517  ind[1] = i + 1 +j*(nx+1);
1518  ind[2] = i + 1 + (j+1)*(nx+1);
1519  ind[3] = i + (j+1)*(nx+1);
1520  elements[k] = new Quadrilateral(ind);
1521  k++;
1522  }
1523  }
1524 
1525  // Sets boundary elements and the corresponding indices of vertices
1526  int m = (nx+1)*ny;
1527  for (i = 0; i < nx; i++)
1528  {
1529  boundary[i] = new Segment(i, i+1, 1);
1530  boundary[nx+i] = new Segment(m+i+1, m+i, 3);
1531  }
1532  m = nx+1;
1533  for (j = 0; j < ny; j++)
1534  {
1535  boundary[2*nx+j] = new Segment((j+1)*m, j*m, 4);
1536  boundary[2*nx+ny+j] = new Segment(j*m+nx, (j+1)*m+nx, 2);
1537  }
1538  }
1539 
1540  // Creates triangular mesh
1541  if (type == Element::TRIANGLE)
1542  {
1543  meshgen = 1;
1544  NumOfVertices = (nx+1) * (ny+1);
1545  NumOfElements = 2 * nx * ny;
1546  NumOfBdrElements = 2 * nx + 2 * ny;
1547  vertices.SetSize(NumOfVertices);
1548  elements.SetSize(NumOfElements);
1549  boundary.SetSize(NumOfBdrElements);
1550  double cx, cy;
1551  int ind[3];
1552 
1553  // Sets vertices and the corresponding coordinates
1554  k = 0;
1555  for (j = 0; j < ny+1; j++)
1556  {
1557  cy = ((double) j / ny) * sy;
1558  for (i = 0; i < nx+1; i++)
1559  {
1560  cx = ((double) i / nx) * sx;
1561  vertices[k](0) = cx;
1562  vertices[k](1) = cy;
1563  k++;
1564  }
1565  }
1566 
1567  // Sets the elements and the corresponding indices of vertices
1568  k = 0;
1569  for (j = 0; j < ny; j++)
1570  {
1571  for (i = 0; i < nx; i++)
1572  {
1573  ind[0] = i + j*(nx+1);
1574  ind[1] = i + 1 + (j+1)*(nx+1);
1575  ind[2] = i + (j+1)*(nx+1);
1576  elements[k] = new Triangle(ind);
1577  k++;
1578  ind[1] = i + 1 + j*(nx+1);
1579  ind[2] = i + 1 + (j+1)*(nx+1);
1580  elements[k] = new Triangle(ind);
1581  k++;
1582  }
1583  }
1584 
1585  // Sets boundary elements and the corresponding indices of vertices
1586  int m = (nx+1)*ny;
1587  for (i = 0; i < nx; i++)
1588  {
1589  boundary[i] = new Segment(i, i+1, 1);
1590  boundary[nx+i] = new Segment(m+i+1, m+i, 3);
1591  }
1592  m = nx+1;
1593  for (j = 0; j < ny; j++)
1594  {
1595  boundary[2*nx+j] = new Segment((j+1)*m, j*m, 4);
1596  boundary[2*nx+ny+j] = new Segment(j*m+nx, (j+1)*m+nx, 2);
1597  }
1598 
1599  MarkTriMeshForRefinement();
1600  }
1601 
1602  CheckElementOrientation();
1603 
1604  if (generate_edges == 1)
1605  {
1606  el_to_edge = new Table;
1607  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1608  GenerateFaces();
1609  CheckBdrElementOrientation();
1610  }
1611  else
1612  NumOfEdges = 0;
1613 
1614  NumOfFaces = 0;
1615 
1616  attributes.Append(1);
1617  bdr_attributes.Append(1); bdr_attributes.Append(2);
1618  bdr_attributes.Append(3); bdr_attributes.Append(4);
1619 }
1620 
1621 void Mesh::Make1D(int n, double sx)
1622 {
1623  int j, ind[1];
1624 
1625  Dim = 1;
1626  spaceDim = 1;
1627 
1628  Init();
1629  InitTables();
1630 
1631  meshgen = 1;
1632 
1633  NumOfVertices = n + 1;
1634  NumOfElements = n;
1635  NumOfBdrElements = 2;
1636  vertices.SetSize(NumOfVertices);
1637  elements.SetSize(NumOfElements);
1638  boundary.SetSize(NumOfBdrElements);
1639 
1640  // Sets vertices and the corresponding coordinates
1641  for (j = 0; j < n+1; j++)
1642  vertices[j](0) = ((double) j / n) * sx;
1643 
1644  // Sets elements and the corresponding indices of vertices
1645  for (j = 0; j < n; j++)
1646  elements[j] = new Segment(j, j+1, 1);
1647 
1648  // Sets the boundary elements
1649  ind[0] = 0;
1650  boundary[0] = new Point(ind, 1);
1651  ind[0] = n;
1652  boundary[1] = new Point(ind, 2);
1653 
1654  NumOfEdges = 0;
1655  NumOfFaces = 0;
1656 
1657  GenerateFaces();
1658 
1659  attributes.Append(1);
1660  bdr_attributes.Append(1); bdr_attributes.Append(2);
1661 }
1662 
1663 Mesh::Mesh(std::istream &input, int generate_edges, int refine, bool fix_orientation)
1664 {
1665  Init();
1666  InitTables();
1667  Load(input, generate_edges, refine, fix_orientation);
1668 }
1669 
1671 {
1672  switch (geom)
1673  {
1674  case Geometry::POINT: return (new Point);
1675  case Geometry::SEGMENT: return (new Segment);
1676  case Geometry::TRIANGLE: return (new Triangle);
1677  case Geometry::SQUARE: return (new Quadrilateral);
1678  case Geometry::CUBE: return (new Hexahedron);
1679  case Geometry::TETRAHEDRON:
1680 #ifdef MFEM_USE_MEMALLOC
1681  return TetMemory.Alloc();
1682 #else
1683  return (new Tetrahedron);
1684 #endif
1685  }
1686 
1687  return NULL;
1688 }
1689 
1691 {
1692  int geom, nv, *v;
1693  Element *el;
1694 
1695  input >> geom;
1696  el = NewElement(geom);
1697  nv = el->GetNVertices();
1698  v = el->GetVertices();
1699  for (int i = 0; i < nv; i++)
1700  input >> v[i];
1701 
1702  return el;
1703 }
1704 
1705 void Mesh::PrintElementWithoutAttr(const Element *el, std::ostream &out)
1706 {
1707  out << el->GetGeometryType();
1708  const int nv = el->GetNVertices();
1709  const int *v = el->GetVertices();
1710  for (int j = 0; j < nv; j++)
1711  out << ' ' << v[j];
1712  out << '\n';
1713 }
1714 
1715 Element *Mesh::ReadElement(std::istream &input)
1716 {
1717  int attr;
1718  Element *el;
1719 
1720  input >> attr;
1721  el = ReadElementWithoutAttr(input);
1722  el->SetAttribute(attr);
1723 
1724  return el;
1725 }
1726 
1727 void Mesh::PrintElement(const Element *el, std::ostream &out)
1728 {
1729  out << el->GetAttribute() << ' ';
1730  PrintElementWithoutAttr(el, out);
1731 }
1732 
1734 {
1735  meshgen = 0;
1736  for (int i = 0; i < NumOfElements; i++)
1737  {
1738  switch (elements[i]->GetType())
1739  {
1740  case Element::SEGMENT:
1741  case Element::TRIANGLE:
1742  case Element::TETRAHEDRON:
1743  meshgen |= 1; break;
1744 
1746  case Element::HEXAHEDRON:
1747  meshgen |= 2;
1748  }
1749  }
1750 }
1751 
1752 // see Tetrahedron::edges
1753 static const int vtk_quadratic_tet[10] =
1754 { 0, 1, 2, 3, 4, 7, 5, 6, 8, 9 };
1755 
1756 // see Hexahedron::edges & Mesh::GenerateFaces
1757 static const int vtk_quadratic_hex[27] =
1758 { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
1759  24, 22, 21, 23, 20, 25, 26 };
1760 
1761 void skip_comment_lines(std::istream &is, const char comment_char)
1762 {
1763  while(1)
1764  {
1765  is >> ws;
1766  if (is.peek() != comment_char)
1767  break;
1768  is.ignore(numeric_limits<streamsize>::max(), '\n');
1769  }
1770 }
1771 
1772 void Mesh::Load(std::istream &input, int generate_edges, int refine,
1773  bool fix_orientation)
1774 {
1775  int i, j, ints[32], n, attr, curved = 0, read_gf = 1;
1776  const int buflen = 1024;
1777  char buf[buflen];
1778 
1779  if (!input)
1780  MFEM_ABORT("Input stream is not open");
1781 
1782  if (NumOfVertices != -1)
1783  {
1784  // Delete the elements.
1785  for (i = 0; i < NumOfElements; i++)
1786  // delete elements[i];
1787  FreeElement(elements[i]);
1788  elements.DeleteAll();
1789 
1790  // Delete the vertices.
1791  vertices.DeleteAll();
1792 
1793  // Delete the boundary elements.
1794  for (i = 0; i < NumOfBdrElements; i++)
1795  // delete boundary[i];
1796  FreeElement(boundary[i]);
1797  boundary.DeleteAll();
1798 
1799  // Delete interior faces (if generated)
1800  for (i = 0; i < faces.Size(); i++)
1801  FreeElement(faces[i]);
1802  faces.DeleteAll();
1803 
1804  faces_info.DeleteAll();
1805 
1806  // Delete the edges (if generated).
1807  DeleteTables();
1808  be_to_edge.DeleteAll();
1809  be_to_face.DeleteAll();
1810 
1811  }
1812 
1813  InitTables();
1814  if (own_nodes) delete Nodes;
1815  Nodes = NULL;
1816  spaceDim = 0;
1817 
1818  string mesh_type;
1819  input >> ws;
1820  getline(input, mesh_type);
1821 
1822  if (mesh_type == "MFEM mesh v1.0")
1823  {
1824  // Read MFEM mesh v1.0 format
1825  string ident;
1826 
1827  // read lines begining with '#' (comments)
1828  skip_comment_lines(input, '#');
1829 
1830  input >> ident; // 'dimension'
1831  input >> Dim;
1832 
1833  skip_comment_lines(input, '#');
1834 
1835  input >> ident; // 'elements'
1836  input >> NumOfElements;
1837  elements.SetSize(NumOfElements);
1838  for (j = 0; j < NumOfElements; j++)
1839  elements[j] = ReadElement(input);
1840 
1841  skip_comment_lines(input, '#');
1842 
1843  input >> ident; // 'boundary'
1844  input >> NumOfBdrElements;
1845  boundary.SetSize(NumOfBdrElements);
1846  for (j = 0; j < NumOfBdrElements; j++)
1847  boundary[j] = ReadElement(input);
1848 
1849  skip_comment_lines(input, '#');
1850 
1851  input >> ident; // 'vertices'
1852  input >> NumOfVertices;
1853  vertices.SetSize(NumOfVertices);
1854 
1855  input >> ws >> ident;
1856  if (ident != "nodes")
1857  {
1858  // read the vertices
1859  spaceDim = atoi(ident.c_str());
1860  for (j = 0; j < NumOfVertices; j++)
1861  for (i = 0; i < spaceDim; i++)
1862  input >> vertices[j](i);
1863  }
1864  else
1865  {
1866  // prepare to read the nodes
1867  input >> ws;
1868  curved = 1;
1869  }
1870  }
1871  else if (mesh_type == "linemesh")
1872  {
1873  int j,p1,p2,a;
1874 
1875  Dim = 1;
1876 
1877  input >> NumOfVertices;
1878  vertices.SetSize(NumOfVertices);
1879  // Sets vertices and the corresponding coordinates
1880  for (j = 0; j < NumOfVertices; j++)
1881  input >> vertices[j](0);
1882 
1883  input >> NumOfElements;
1884  elements.SetSize(NumOfElements);
1885  // Sets elements and the corresponding indices of vertices
1886  for (j = 0; j < NumOfElements; j++)
1887  {
1888  input >> a >> p1 >> p2;
1889  elements[j] = new Segment(p1-1, p2-1, a);
1890  }
1891 
1892  int ind[1];
1893  input >> NumOfBdrElements;
1894  boundary.SetSize(NumOfBdrElements);
1895  for (j = 0; j < NumOfBdrElements; j++)
1896  {
1897  input >> a >> ind[0];
1898  ind[0]--;
1899  boundary[j] = new Point(ind,a);
1900  }
1901  }
1902  else if (mesh_type == "areamesh2" || mesh_type == "curved_areamesh2")
1903  {
1904  // Read planar mesh in Netgen format.
1905  Dim = 2;
1906 
1907  if (mesh_type == "curved_areamesh2")
1908  curved = 1;
1909 
1910  // Read the boundary elements.
1911  input >> NumOfBdrElements;
1912  boundary.SetSize(NumOfBdrElements);
1913  for (i = 0; i < NumOfBdrElements; i++)
1914  {
1915  input >> attr
1916  >> ints[0] >> ints[1];
1917  ints[0]--; ints[1]--;
1918  boundary[i] = new Segment(ints, attr);
1919  }
1920 
1921  // Read the elements.
1922  input >> NumOfElements;
1923  elements.SetSize(NumOfElements);
1924  for (i = 0; i < NumOfElements; i++)
1925  {
1926  input >> attr >> n;
1927  for (j = 0; j < n; j++)
1928  {
1929  input >> ints[j];
1930  ints[j]--;
1931  }
1932  switch (n)
1933  {
1934  case 2:
1935  elements[i] = new Segment(ints, attr);
1936  break;
1937  case 3:
1938  elements[i] = new Triangle(ints, attr);
1939  break;
1940  case 4:
1941  elements[i] = new Quadrilateral(ints, attr);
1942  break;
1943  }
1944  }
1945 
1946  if (!curved)
1947  {
1948  // Read the vertices.
1949  input >> NumOfVertices;
1950  vertices.SetSize(NumOfVertices);
1951  for (i = 0; i < NumOfVertices; i++)
1952  for (j = 0; j < Dim; j++)
1953  input >> vertices[i](j);
1954  }
1955  else
1956  {
1957  input >> NumOfVertices;
1958  vertices.SetSize(NumOfVertices);
1959  input >> ws;
1960  }
1961  }
1962  else if (mesh_type == "NETGEN" || mesh_type == "NETGEN_Neutral_Format")
1963  {
1964  // Read a netgen format mesh of tetrahedra.
1965  Dim = 3;
1966 
1967  // Read the vertices
1968  input >> NumOfVertices;
1969 
1970  vertices.SetSize(NumOfVertices);
1971  for (i = 0; i < NumOfVertices; i++)
1972  for (j = 0; j < Dim; j++)
1973  input >> vertices[i](j);
1974 
1975  // Read the elements
1976  input >> NumOfElements;
1977  elements.SetSize(NumOfElements);
1978  for (i = 0; i < NumOfElements; i++)
1979  {
1980  input >> attr;
1981  for (j = 0; j < 4; j++)
1982  {
1983  input >> ints[j];
1984  ints[j]--;
1985  }
1986 #ifdef MFEM_USE_MEMALLOC
1987  Tetrahedron *tet;
1988  tet = TetMemory.Alloc();
1989  tet->SetVertices(ints);
1990  tet->SetAttribute(attr);
1991  elements[i] = tet;
1992 #else
1993  elements[i] = new Tetrahedron(ints, attr);
1994 #endif
1995  }
1996 
1997  // Read the boundary information.
1998  input >> NumOfBdrElements;
1999  boundary.SetSize(NumOfBdrElements);
2000  for (i = 0; i < NumOfBdrElements; i++)
2001  {
2002  input >> attr;
2003  for (j = 0; j < 3; j++)
2004  {
2005  input >> ints[j];
2006  ints[j]--;
2007  }
2008  boundary[i] = new Triangle(ints, attr);
2009  }
2010  }
2011  else if (mesh_type == "TrueGrid")
2012  {
2013  // Reading TrueGrid mesh.
2014 
2015  // TODO: find the actual dimension
2016  Dim = 3;
2017 
2018  if (Dim == 2)
2019  {
2020  int vari;
2021  double varf;
2022 
2023  input >> vari >> NumOfVertices >> vari >> vari >> NumOfElements;
2024  input.getline(buf, buflen);
2025  input.getline(buf, buflen);
2026  input >> vari;
2027  input.getline(buf, buflen);
2028  input.getline(buf, buflen);
2029  input.getline(buf, buflen);
2030 
2031  // Read the vertices.
2032  vertices.SetSize(NumOfVertices);
2033  for (i = 0; i < NumOfVertices; i++)
2034  {
2035  input >> vari >> varf >> vertices[i](0) >> vertices[i](1);
2036  input.getline(buf, buflen);
2037  }
2038 
2039  // Read the elements.
2040  elements.SetSize(NumOfElements);
2041  for (i = 0; i < NumOfElements; i++)
2042  {
2043  input >> vari >> attr;
2044  for (j = 0; j < 4; j++)
2045  {
2046  input >> ints[j];
2047  ints[j]--;
2048  }
2049  input.getline(buf, buflen);
2050  input.getline(buf, buflen);
2051  elements[i] = new Quadrilateral(ints, attr);
2052  }
2053  }
2054  else if (Dim == 3)
2055  {
2056  int vari;
2057  double varf;
2058  input >> vari >> NumOfVertices >> NumOfElements;
2059  input.getline(buf, buflen);
2060  input.getline(buf, buflen);
2061  input >> vari >> vari >> NumOfBdrElements;
2062  input.getline(buf, buflen);
2063  input.getline(buf, buflen);
2064  input.getline(buf, buflen);
2065  // Read the vertices.
2066  vertices.SetSize(NumOfVertices);
2067  for (i = 0; i < NumOfVertices; i++)
2068  {
2069  input >> vari >> varf >> vertices[i](0) >> vertices[i](1)
2070  >> vertices[i](2);
2071  input.getline(buf, buflen);
2072  }
2073  // Read the elements.
2074  elements.SetSize(NumOfElements);
2075  for (i = 0; i < NumOfElements; i++)
2076  {
2077  input >> vari >> attr;
2078  for (j = 0; j < 8; j++)
2079  {
2080  input >> ints[j];
2081  ints[j]--;
2082  }
2083  input.getline(buf, buflen);
2084  elements[i] = new Hexahedron(ints, attr);
2085  }
2086  // Read the boundary elements.
2087  boundary.SetSize(NumOfBdrElements);
2088  for (i = 0; i < NumOfBdrElements; i++)
2089  {
2090  input >> attr;
2091  for (j = 0; j < 4; j++)
2092  {
2093  input >> ints[j];
2094  ints[j]--;
2095  }
2096  input.getline(buf, buflen);
2097  boundary[i] = new Quadrilateral(ints, attr);
2098  }
2099  }
2100  }
2101  else if (mesh_type == "# vtk DataFile Version 3.0" ||
2102  mesh_type == "# vtk DataFile Version 2.0")
2103  {
2104  // Reading VTK mesh
2105 
2106  string buff;
2107  getline(input, buff); // comment line
2108  getline(input, buff);
2109  if (buff != "ASCII")
2110  {
2111  mfem_error("Mesh::Load : VTK mesh is not in ASCII format!");
2112  return;
2113  }
2114  getline(input, buff);
2115  if (buff != "DATASET UNSTRUCTURED_GRID")
2116  {
2117  mfem_error("Mesh::Load : VTK mesh is not UNSTRUCTURED_GRID!");
2118  return;
2119  }
2120 
2121  // Read the points, skipping optional sections such as the FIELD data from
2122  // VisIt's VTK export (or from Mesh::PrintVTK with field_data==1).
2123  do
2124  {
2125  input >> buff;
2126  if (!input.good())
2127  mfem_error("Mesh::Load : VTK mesh does not have POINTS data!");
2128  }
2129  while (buff != "POINTS");
2130  int np = 0;
2131  Vector points;
2132  {
2133  input >> np >> ws;
2134  points.SetSize(3*np);
2135  getline(input, buff); // "double"
2136  for (i = 0; i < points.Size(); i++)
2137  input >> points(i);
2138  }
2139 
2140  // Read the cells
2141  NumOfElements = n = 0;
2142  Array<int> cells_data;
2143  input >> ws >> buff;
2144  if (buff == "CELLS")
2145  {
2146  input >> NumOfElements >> n >> ws;
2147  cells_data.SetSize(n);
2148  for (i = 0; i < n; i++)
2149  input >> cells_data[i];
2150  }
2151 
2152  // Read the cell types
2153  Dim = 0;
2154  int order = 1;
2155  input >> ws >> buff;
2156  if (buff == "CELL_TYPES")
2157  {
2158  input >> NumOfElements;
2159  elements.SetSize(NumOfElements);
2160  for (j = i = 0; i < NumOfElements; i++)
2161  {
2162  int ct;
2163  input >> ct;
2164  switch (ct)
2165  {
2166  case 5: // triangle
2167  Dim = 2;
2168  elements[i] = new Triangle(&cells_data[j+1]);
2169  break;
2170  case 9: // quadrilateral
2171  Dim = 2;
2172  elements[i] = new Quadrilateral(&cells_data[j+1]);
2173  break;
2174  case 10: // tetrahedron
2175  Dim = 3;
2176 #ifdef MFEM_USE_MEMALLOC
2177  elements[i] = TetMemory.Alloc();
2178  elements[i]->SetVertices(&cells_data[j+1]);
2179 #else
2180  elements[i] = new Tetrahedron(&cells_data[j+1]);
2181 #endif
2182  break;
2183  case 12: // hexahedron
2184  Dim = 3;
2185  elements[i] = new Hexahedron(&cells_data[j+1]);
2186  break;
2187 
2188  case 22: // quadratic triangle
2189  Dim = 2;
2190  order = 2;
2191  elements[i] = new Triangle(&cells_data[j+1]);
2192  break;
2193  case 28: // biquadratic quadrilateral
2194  Dim = 2;
2195  order = 2;
2196  elements[i] = new Quadrilateral(&cells_data[j+1]);
2197  break;
2198  case 24: // quadratic tetrahedron
2199  Dim = 3;
2200  order = 2;
2201 #ifdef MFEM_USE_MEMALLOC
2202  elements[i] = TetMemory.Alloc();
2203  elements[i]->SetVertices(&cells_data[j+1]);
2204 #else
2205  elements[i] = new Tetrahedron(&cells_data[j+1]);
2206 #endif
2207  break;
2208  case 29: // triquadratic hexahedron
2209  Dim = 3;
2210  order = 2;
2211  elements[i] = new Hexahedron(&cells_data[j+1]);
2212  break;
2213  default:
2214  cerr << "Mesh::Load : VTK mesh : cell type " << ct
2215  << " is not supported!" << endl;
2216  mfem_error();
2217  return;
2218  }
2219  j += cells_data[j] + 1;
2220  }
2221  }
2222 
2223  // Read attributes
2224  streampos sp = input.tellg();
2225  input >> ws >> buff;
2226  if (buff == "CELL_DATA")
2227  {
2228  input >> n >> ws;
2229  getline(input, buff);
2230  if (buff == "SCALARS material int" || buff == "SCALARS material float")
2231  {
2232  getline(input, buff); // "LOOKUP_TABLE default"
2233  for (i = 0; i < NumOfElements; i++)
2234  {
2235  input >> attr;
2236  elements[i]->SetAttribute(attr);
2237  }
2238  }
2239  else
2240  input.seekg(sp);
2241  }
2242  else
2243  input.seekg(sp);
2244 
2245  if (order == 1)
2246  {
2247  cells_data.DeleteAll();
2248  NumOfVertices = np;
2249  vertices.SetSize(np);
2250  for (i = 0; i < np; i++)
2251  {
2252  vertices[i](0) = points(3*i+0);
2253  vertices[i](1) = points(3*i+1);
2254  vertices[i](2) = points(3*i+2);
2255  }
2256  points.Destroy();
2257 
2258  // No boundary is defined in a VTK mesh
2259  NumOfBdrElements = 0;
2260  }
2261  else if (order == 2)
2262  {
2263  curved = 1;
2264 
2265  // generate new enumeration for the vertices
2266  Array<int> pts_dof(np);
2267  pts_dof = -1;
2268  for (n = i = 0; i < NumOfElements; i++)
2269  {
2270  int *v = elements[i]->GetVertices();
2271  int nv = elements[i]->GetNVertices();
2272  for (j = 0; j < nv; j++)
2273  if (pts_dof[v[j]] == -1)
2274  pts_dof[v[j]] = n++;
2275  }
2276  // keep the original ordering of the vertices
2277  for (n = i = 0; i < np; i++)
2278  if (pts_dof[i] != -1)
2279  pts_dof[i] = n++;
2280  // update the element vertices
2281  for (i = 0; i < NumOfElements; i++)
2282  {
2283  int *v = elements[i]->GetVertices();
2284  int nv = elements[i]->GetNVertices();
2285  for (j = 0; j < nv; j++)
2286  v[j] = pts_dof[v[j]];
2287  }
2288  // Define the 'vertices' from the 'points' through the 'pts_dof' map
2289  NumOfVertices = n;
2290  vertices.SetSize(n);
2291  for (i = 0; i < np; i++)
2292  {
2293  if ((j = pts_dof[i]) != -1)
2294  {
2295  vertices[j](0) = points(3*i+0);
2296  vertices[j](1) = points(3*i+1);
2297  vertices[j](2) = points(3*i+2);
2298  }
2299  }
2300 
2301  // No boundary is defined in a VTK mesh
2302  NumOfBdrElements = 0;
2303 
2304  // Generate faces and edges so that we can define quadratic
2305  // FE space on the mesh
2306 
2307  // Generate faces
2308  if (Dim > 2)
2309  {
2310  GetElementToFaceTable();
2311  GenerateFaces();
2312  }
2313  else
2314  NumOfFaces = 0;
2315 
2316  // Generate edges
2317  el_to_edge = new Table;
2318  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2319  if (Dim == 2)
2320  GenerateFaces(); // 'Faces' in 2D refers to the edges
2321 
2322  // Define quadratic FE space
2324  FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim);
2325  Nodes = new GridFunction(fes);
2326  Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
2327  own_nodes = 1;
2328 
2329  // Map vtk points to edge/face/element dofs
2330  Array<int> dofs;
2331  for (n = i = 0; i < NumOfElements; i++)
2332  {
2333  fes->GetElementDofs(i, dofs);
2334  const int *vtk_mfem;
2335  switch (elements[i]->GetGeometryType())
2336  {
2337  case Geometry::TRIANGLE:
2338  case Geometry::SQUARE:
2339  vtk_mfem = vtk_quadratic_hex; break; // identity map
2340  case Geometry::TETRAHEDRON:
2341  vtk_mfem = vtk_quadratic_tet; break;
2342  case Geometry::CUBE:
2343  default:
2344  vtk_mfem = vtk_quadratic_hex; break;
2345  }
2346 
2347  for (n++, j = 0; j < dofs.Size(); j++, n++)
2348  {
2349  if (pts_dof[cells_data[n]] == -1)
2350  {
2351  pts_dof[cells_data[n]] = dofs[vtk_mfem[j]];
2352  }
2353  else
2354  {
2355  if (pts_dof[cells_data[n]] != dofs[vtk_mfem[j]])
2356  mfem_error("Mesh::Load : VTK mesh : "
2357  "inconsistent quadratic mesh!");
2358  }
2359  }
2360  }
2361 
2362  // Define the 'Nodes' from the 'points' through the 'pts_dof' map
2363  for (i = 0; i < np; i++)
2364  {
2365  dofs.SetSize(1);
2366  if ((dofs[0] = pts_dof[i]) != -1)
2367  {
2368  fes->DofsToVDofs(dofs);
2369  for (j = 0; j < dofs.Size(); j++)
2370  (*Nodes)(dofs[j]) = points(3*i+j);
2371  }
2372  }
2373 
2374  read_gf = 0;
2375  }
2376  }
2377  else if (mesh_type == "MFEM NURBS mesh v1.0")
2378  {
2379  NURBSext = new NURBSExtension(input);
2380 
2381  Dim = NURBSext->Dimension();
2382  NumOfVertices = NURBSext->GetNV();
2383  NumOfElements = NURBSext->GetNE();
2384  NumOfBdrElements = NURBSext->GetNBE();
2385 
2386  NURBSext->GetElementTopo(elements);
2387  NURBSext->GetBdrElementTopo(boundary);
2388 
2389  vertices.SetSize(NumOfVertices);
2390  curved = 1;
2391  if (NURBSext->HavePatches())
2392  {
2393  NURBSFECollection *fec = new NURBSFECollection(NURBSext->GetOrder());
2394  FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim,
2396  Nodes = new GridFunction(fes);
2397  Nodes->MakeOwner(fec);
2398  NURBSext->SetCoordsFromPatches(*Nodes);
2399  own_nodes = 1;
2400  read_gf = 0;
2401  int vd = Nodes->VectorDim();
2402  for (i = 0; i < vd; i++)
2403  {
2404  Vector vert_val;
2405  Nodes->GetNodalValues(vert_val, i+1);
2406  for (j = 0; j < NumOfVertices; j++)
2407  vertices[j](i) = vert_val(j);
2408  }
2409  }
2410  else
2411  read_gf = 1;
2412  }
2413  else if (mesh_type == "MFEM INLINE mesh v1.0")
2414  {
2415  // Initialize to negative numbers so that we know if they've
2416  // been set. We're using Element::POINT as our flag, since
2417  // we're not going to make a 0D mesh, ever.
2418  int nx = -1;
2419  int ny = -1;
2420  int nz = -1;
2421  double sx = -1.0;
2422  double sy = -1.0;
2423  double sz = -1.0;
2425 
2426  while (true)
2427  {
2428  skip_comment_lines(input, '#');
2429  // Break out if we reached the end of the file after
2430  // gobbling up the whitespace and comments after the last keyword.
2431  if (!input.good())
2432  {
2433  break;
2434  }
2435 
2436  // Read the next keyword
2437  std::string name;
2438  input >> name;
2439  input >> std::ws;
2440  // Make sure there's an equal sign
2441  if (input.get() != '=')
2442  {
2443  std::ostringstream os;
2444  os << "Mesh::Load : Inline mesh expected '=' after keyword " << name;
2445  mfem_error(os.str().c_str());
2446  }
2447  input >> std::ws;
2448 
2449  if (name == "nx")
2450  {
2451  input >> nx;
2452  }
2453  else if (name == "ny")
2454  {
2455  input >> ny;
2456  }
2457  else if (name == "nz")
2458  {
2459  input >> nz;
2460  }
2461  else if (name == "sx")
2462  {
2463  input >> sx;
2464  }
2465  else if (name == "sy")
2466  {
2467  input >> sy;
2468  }
2469  else if (name == "sz")
2470  {
2471  input >> sz;
2472  }
2473  else if (name == "type")
2474  {
2475  std::string eltype;
2476  input >> eltype;
2477  if (eltype == "segment")
2478  {
2479  type = Element::SEGMENT;
2480  }
2481  else if (eltype == "quad")
2482  {
2483  type = Element::QUADRILATERAL;
2484  }
2485  else if (eltype == "tri")
2486  {
2487  type = Element::TRIANGLE;
2488  }
2489  else if (eltype == "hex")
2490  {
2491  type = Element::HEXAHEDRON;
2492  }
2493  else if (eltype == "tet")
2494  {
2495  type = Element::TETRAHEDRON;
2496  }
2497  else
2498  {
2499  std::ostringstream os;
2500  os << "Mesh::Load : unrecognized element type (read '" << eltype
2501  << "') in inline mesh format. Allowed: segment, tri, tet, quad, hex";
2502  mfem_error(os.str().c_str());
2503  }
2504  }
2505  else
2506  {
2507  std::ostringstream os;
2508  os << "Mesh::Load : unrecognized keyword (" << name
2509  << ") in inline mesh format. Allowed: nx, ny, nz, type, sx, sy, sz";
2510  mfem_error(os.str().c_str());
2511  }
2512 
2513  input >> std::ws;
2514  // Allow an optional semi-colon at the end of each line.
2515  if (input.peek() == ';')
2516  {
2517  input.get();
2518  }
2519 
2520  // Done reading file
2521  if (!input)
2522  {
2523  break;
2524  }
2525  }
2526 
2527  // Now make the mesh.
2528  if (type == Element::SEGMENT)
2529  {
2530  if (nx < 0 || sx < 0.0)
2531  {
2532  std::ostringstream os;
2533  os << "Mesh::Load : invalid 1D inline mesh format, all values must be "
2534  "positive\n"
2535  << " nx = " << nx << "\n"
2536  << " sx = " << sx << "\n";
2537  mfem_error(os.str().c_str());
2538  }
2539  Make1D(nx, sx);
2540  }
2541  else if (type == Element::TRIANGLE || type == Element::QUADRILATERAL)
2542  {
2543  if (nx < 0 || ny < 0 || sx < 0.0 || sy < 0.0)
2544  {
2545  std::ostringstream os;
2546  os << "Mesh::Load : invalid 2D inline mesh format, all values must be "
2547  "positive\n"
2548  << " nx = " << nx << "\n"
2549  << " ny = " << ny << "\n"
2550  << " sx = " << sx << "\n"
2551  << " sy = " << sy << "\n";
2552  mfem_error(os.str().c_str());
2553  }
2554  Make2D(nx, ny, type, generate_edges, sx, sy);
2555  }
2556  else if (type == Element::TETRAHEDRON || type == Element::HEXAHEDRON)
2557  {
2558  if (nx < 0 || ny < 0 || nz < 0 || sx < 0.0 || sy < 0.0 || sz < 0.0)
2559  {
2560  std::ostringstream os;
2561  os << "Mesh::Load : invalid 3D inline mesh format, all values must be "
2562  "positive\n"
2563  << " nx = " << nx << "\n"
2564  << " ny = " << ny << "\n"
2565  << " nz = " << nz << "\n"
2566  << " sx = " << sx << "\n"
2567  << " sy = " << sy << "\n"
2568  << " sz = " << sz << "\n";
2569  mfem_error(os.str().c_str());
2570  }
2571  Make3D(nx, ny, nz, type, generate_edges, sx, sy, sz);
2572  }
2573  else
2574  {
2575  mfem_error("Mesh::Load : For inline mesh, must specify an "
2576  "element type = [segment, tri, quad, tet, hex]");
2577  }
2578  return; // done with inline mesh construction
2579  }
2580  else
2581  {
2582  MFEM_ABORT("Unknown input mesh format");
2583  return;
2584  }
2585 
2586  // at this point the following should be defined:
2587  // 1) Dim
2588  // 2) NumOfElements, elements
2589  // 3) NumOfBdrElements, boundary
2590  // 4) NumOfVertices, with allocated space in vertices
2591  // 5) curved
2592  // 5a) if curved == 0, vertices must be defined
2593  // 5b) if curved != 0 and read_gf != 0,
2594  // 'input' must point to a GridFunction
2595  // 5c) if curved != 0 and read_gf == 0,
2596  // vertices and Nodes must be defined
2597 
2598  if (spaceDim == 0)
2599  spaceDim = Dim;
2600 
2601  // set the mesh type ('meshgen')
2602  SetMeshGen();
2603 
2604  if (NumOfBdrElements == 0 && Dim > 2)
2605  {
2606  // in 3D, generate boundary elements before we 'MarkForRefinement'
2607  GetElementToFaceTable();
2608  GenerateFaces();
2609  GenerateBoundaryElements();
2610  }
2611 
2612  if (!curved)
2613  {
2614  // check and fix element orientation
2615  CheckElementOrientation(fix_orientation);
2616 
2617  if (refine)
2618  MarkForRefinement();
2619  }
2620 
2621  if (Dim == 1)
2622  GenerateFaces();
2623 
2624  // generate the faces
2625  if (Dim > 2)
2626  {
2627  GetElementToFaceTable();
2628  GenerateFaces();
2629  // check and fix boundary element orientation
2630  if ( !(curved && (meshgen & 1)) )
2631  CheckBdrElementOrientation();
2632  }
2633  else
2634  NumOfFaces = 0;
2635 
2636  // generate edges if requested
2637  if (Dim > 1 && generate_edges == 1)
2638  {
2639  el_to_edge = new Table;
2640  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2641  if (Dim == 2)
2642  {
2643  GenerateFaces(); // 'Faces' in 2D refers to the edges
2644  if (NumOfBdrElements == 0)
2645  GenerateBoundaryElements();
2646  // check and fix boundary element orientation
2647  if ( !(curved && (meshgen & 1)) )
2648  CheckBdrElementOrientation();
2649  }
2650  c_el_to_edge = NULL;
2651  }
2652  else
2653  NumOfEdges = 0;
2654 
2655  // generate the arrays 'attributes' and ' bdr_attributes'
2656  SetAttributes();
2657 
2658  if (curved)
2659  {
2660  if (read_gf)
2661  {
2662  Nodes = new GridFunction(this, input);
2663  own_nodes = 1;
2664  spaceDim = Nodes->VectorDim();
2665  // Set the 'vertices' from the 'Nodes'
2666  for (i = 0; i < spaceDim; i++)
2667  {
2668  Vector vert_val;
2669  Nodes->GetNodalValues(vert_val, i+1);
2670  for (j = 0; j < NumOfVertices; j++)
2671  vertices[j](i) = vert_val(j);
2672  }
2673  }
2674 
2675  // Check orientation and mark edges; only for triangles / tets
2676  if (meshgen & 1)
2677  {
2678  DSTable *old_v_to_v = NULL;
2679  Table *old_elem_vert = NULL;
2680  if (fix_orientation || refine)
2681  {
2682  PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
2683  }
2684 
2685  // check orientation and mark for refinement using just vertices
2686  // (i.e. higher order curvature is not used)
2687  CheckElementOrientation(fix_orientation);
2688  if (refine)
2689  MarkForRefinement(); // changes topology!
2690 
2691  if (fix_orientation || refine)
2692  {
2693  DoNodeReorder(old_v_to_v, old_elem_vert);
2694  delete old_elem_vert;
2695  delete old_v_to_v;
2696  }
2697  }
2698  }
2699 }
2700 
2701 Mesh::Mesh(Mesh *mesh_array[], int num_pieces)
2702 {
2703  int i, j, ie, ib, iv, *v, nv;
2704  Element *el;
2705  Mesh *m;
2706 
2707  Init();
2708  InitTables();
2709 
2710  Dim = mesh_array[0]->Dimension();
2711  spaceDim = mesh_array[0]->SpaceDimension();
2712 
2713  if (mesh_array[0]->NURBSext)
2714  {
2715  // assuming the pieces form a partition of a NURBS mesh
2716  NURBSext = new NURBSExtension(mesh_array, num_pieces);
2717 
2718  NumOfVertices = NURBSext->GetNV();
2719  NumOfElements = NURBSext->GetNE();
2720 
2721  NURBSext->GetElementTopo(elements);
2722 
2723  // NumOfBdrElements = NURBSext->GetNBE();
2724  // NURBSext->GetBdrElementTopo(boundary);
2725 
2726  Array<int> lvert_vert, lelem_elem;
2727 
2728  // Here, for visualization purposes, we copy the boundary elements from
2729  // the individual pieces which include the interior boundaries.
2730  // This creates 'boundary' array that is different from the one generated
2731  // by the NURBSExtension which, in particular, makes the boundary-dof
2732  // table invalid. This, in turn, causes GetBdrElementTransformation to
2733  // not function properly.
2734  NumOfBdrElements = 0;
2735  for (i = 0; i < num_pieces; i++)
2736  NumOfBdrElements += mesh_array[i]->GetNBE();
2737  boundary.SetSize(NumOfBdrElements);
2738  vertices.SetSize(NumOfVertices);
2739  ib = 0;
2740  for (i = 0; i < num_pieces; i++)
2741  {
2742  m = mesh_array[i];
2743  m->NURBSext->GetVertexLocalToGlobal(lvert_vert);
2744  m->NURBSext->GetElementLocalToGlobal(lelem_elem);
2745  // copy the element attributes
2746  for (j = 0; j < m->GetNE(); j++)
2747  elements[lelem_elem[j]]->SetAttribute(m->GetAttribute(j));
2748  // copy the boundary
2749  for (j = 0; j < m->GetNBE(); j++)
2750  {
2751  el = m->GetBdrElement(j)->Duplicate(this);
2752  v = el->GetVertices();
2753  nv = el->GetNVertices();
2754  for (int k = 0; k < nv; k++)
2755  v[k] = lvert_vert[v[k]];
2756  boundary[ib++] = el;
2757  }
2758  // copy the vertices
2759  for (j = 0; j < m->GetNV(); j++)
2760  vertices[lvert_vert[j]].SetCoords(m->GetVertex(j));
2761  }
2762  }
2763  else // not a NURBS mesh
2764  {
2765  NumOfElements = 0;
2766  NumOfBdrElements = 0;
2767  NumOfVertices = 0;
2768  for (i = 0; i < num_pieces; i++)
2769  {
2770  m = mesh_array[i];
2771  NumOfElements += m->GetNE();
2772  NumOfBdrElements += m->GetNBE();
2773  NumOfVertices += m->GetNV();
2774  }
2775  elements.SetSize(NumOfElements);
2776  boundary.SetSize(NumOfBdrElements);
2777  vertices.SetSize(NumOfVertices);
2778  ie = ib = iv = 0;
2779  for (i = 0; i < num_pieces; i++)
2780  {
2781  m = mesh_array[i];
2782  // copy the elements
2783  for (j = 0; j < m->GetNE(); j++)
2784  {
2785  el = m->GetElement(j)->Duplicate(this);
2786  v = el->GetVertices();
2787  nv = el->GetNVertices();
2788  for (int k = 0; k < nv; k++)
2789  v[k] += iv;
2790  elements[ie++] = el;
2791  }
2792  // copy the boundary elements
2793  for (j = 0; j < m->GetNBE(); j++)
2794  {
2795  el = m->GetBdrElement(j)->Duplicate(this);
2796  v = el->GetVertices();
2797  nv = el->GetNVertices();
2798  for (int k = 0; k < nv; k++)
2799  v[k] += iv;
2800  boundary[ib++] = el;
2801  }
2802  // copy the vertices
2803  for (j = 0; j < m->GetNV(); j++)
2804  vertices[iv++].SetCoords(m->GetVertex(j));
2805  }
2806  }
2807 
2808  // set the mesh type ('meshgen')
2809  meshgen = 0;
2810  for (i = 0; i < num_pieces; i++)
2811  meshgen |= mesh_array[i]->MeshGenerator();
2812 
2813  // generate faces
2814  if (Dim > 2)
2815  {
2816  GetElementToFaceTable();
2817  GenerateFaces();
2818  }
2819  else
2820  NumOfFaces = 0;
2821 
2822  // generate edges
2823  if (Dim > 1)
2824  {
2825  el_to_edge = new Table;
2826  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2827  if (Dim == 2)
2828  GenerateFaces(); // 'Faces' in 2D refers to the edges
2829  }
2830  else
2831  NumOfEdges = 0;
2832 
2833  // generate the arrays 'attributes' and ' bdr_attributes'
2834  SetAttributes();
2835 
2836  // copy the nodes (curvilinear meshes)
2837  GridFunction *g = mesh_array[0]->GetNodes();
2838  if (g)
2839  {
2840  Array<GridFunction *> gf_array(num_pieces);
2841  for (i = 0; i < num_pieces; i++)
2842  gf_array[i] = mesh_array[i]->GetNodes();
2843  Nodes = new GridFunction(this, gf_array, num_pieces);
2844  own_nodes = 1;
2845  }
2846 
2847 #ifdef MFEM_DEBUG
2848  CheckElementOrientation(false);
2849  CheckBdrElementOrientation(false);
2850 #endif
2851 }
2852 
2854 {
2855  if (NURBSext == NULL)
2856  mfem_error("Mesh::KnotInsert : Not a NURBS mesh!");
2857 
2858  if (kv.Size() != NURBSext->GetNKV())
2859  mfem_error("Mesh::KnotInsert : KnotVector array size mismatch!");
2860 
2861  NURBSext->ConvertToPatches(*Nodes);
2862 
2863  NURBSext->KnotInsert(kv);
2864 
2865  UpdateNURBS();
2866 }
2867 
2869 {
2870  // do not check for NURBSext since this method is protected
2871  NURBSext->ConvertToPatches(*Nodes);
2872 
2873  NURBSext->UniformRefinement();
2874 
2875  UpdateNURBS();
2876 }
2877 
2879 {
2880  if (NURBSext == NULL)
2881  mfem_error("Mesh::DegreeElevate : Not a NURBS mesh!");
2882 
2883  NURBSext->ConvertToPatches(*Nodes);
2884 
2885  NURBSext->DegreeElevate(t);
2886 
2887  NURBSFECollection *nurbs_fec =
2888  dynamic_cast<NURBSFECollection *>(Nodes->OwnFEC());
2889  if (!nurbs_fec)
2890  mfem_error("Mesh::DegreeElevate");
2891  nurbs_fec->UpdateOrder(nurbs_fec->GetOrder() + t);
2892 
2893  UpdateNURBS();
2894 }
2895 
2897 {
2898  NURBSext->SetKnotsFromPatches();
2899 
2900  Dim = NURBSext->Dimension();
2901  spaceDim = Dim;
2902 
2903  if (NumOfElements != NURBSext->GetNE())
2904  {
2905  for (int i = 0; i < elements.Size(); i++)
2906  FreeElement(elements[i]);
2907  NumOfElements = NURBSext->GetNE();
2908  NURBSext->GetElementTopo(elements);
2909  }
2910 
2911  if (NumOfBdrElements != NURBSext->GetNBE())
2912  {
2913  for (int i = 0; i < boundary.Size(); i++)
2914  FreeElement(boundary[i]);
2915  NumOfBdrElements = NURBSext->GetNBE();
2916  NURBSext->GetBdrElementTopo(boundary);
2917  }
2918 
2919  Nodes->FESpace()->Update();
2920  Nodes->Update();
2921  NURBSext->SetCoordsFromPatches(*Nodes);
2922 
2923  if (NumOfVertices != NURBSext->GetNV())
2924  {
2925  NumOfVertices = NURBSext->GetNV();
2926  vertices.SetSize(NumOfVertices);
2927  int vd = Nodes->VectorDim();
2928  for (int i = 0; i < vd; i++)
2929  {
2930  Vector vert_val;
2931  Nodes->GetNodalValues(vert_val, i+1);
2932  for (int j = 0; j < NumOfVertices; j++)
2933  vertices[j](i) = vert_val(j);
2934  }
2935  }
2936 
2937  if (el_to_edge)
2938  {
2939  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2940  if (Dim == 2)
2941  GenerateFaces();
2942  }
2943 
2944  if (el_to_face)
2945  {
2946  GetElementToFaceTable();
2947  GenerateFaces();
2948  }
2949 }
2950 
2951 void Mesh::LoadPatchTopo(std::istream &input, Array<int> &edge_to_knot)
2952 {
2953  Init();
2954  InitTables();
2955 
2956  int j;
2957 
2958  // Read MFEM NURBS mesh v1.0 format
2959  string ident;
2960 
2961  skip_comment_lines(input, '#');
2962 
2963  input >> ident; // 'dimension'
2964  input >> Dim;
2965  spaceDim = Dim;
2966 
2967  skip_comment_lines(input, '#');
2968 
2969  input >> ident; // 'elements'
2970  input >> NumOfElements;
2971  elements.SetSize(NumOfElements);
2972  for (j = 0; j < NumOfElements; j++)
2973  elements[j] = ReadElement(input);
2974 
2975  skip_comment_lines(input, '#');
2976 
2977  input >> ident; // 'boundary'
2978  input >> NumOfBdrElements;
2979  boundary.SetSize(NumOfBdrElements);
2980  for (j = 0; j < NumOfBdrElements; j++)
2981  boundary[j] = ReadElement(input);
2982 
2983  skip_comment_lines(input, '#');
2984 
2985  input >> ident; // 'edges'
2986  input >> NumOfEdges;
2987  edge_vertex = new Table(NumOfEdges, 2);
2988  edge_to_knot.SetSize(NumOfEdges);
2989  for (j = 0; j < NumOfEdges; j++)
2990  {
2991  int *v = edge_vertex->GetRow(j);
2992  input >> edge_to_knot[j] >> v[0] >> v[1];
2993  if (v[0] > v[1])
2994  edge_to_knot[j] = -1 - edge_to_knot[j];
2995  }
2996 
2997  skip_comment_lines(input, '#');
2998 
2999  input >> ident; // 'vertices'
3000  input >> NumOfVertices;
3001  vertices.SetSize(0);
3002 
3003  meshgen = 2;
3004 
3005  // generate the faces
3006  if (Dim > 2)
3007  {
3008  GetElementToFaceTable();
3009  GenerateFaces();
3010  if (NumOfBdrElements == 0)
3011  GenerateBoundaryElements();
3012  CheckBdrElementOrientation();
3013  }
3014  else
3015  NumOfFaces = 0;
3016 
3017  // generate edges
3018  if (Dim > 1)
3019  {
3020  el_to_edge = new Table;
3021  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3022  if (Dim < 3)
3023  {
3024  GenerateFaces();
3025  if (NumOfBdrElements == 0)
3026  GenerateBoundaryElements();
3027  CheckBdrElementOrientation();
3028  }
3029  }
3030  else
3031  NumOfEdges = 0;
3032 
3033  // generate the arrays 'attributes' and ' bdr_attributes'
3034  SetAttributes();
3035 }
3036 
3038 {
3039  if (p.Size() >= v.Size())
3040  {
3041  for (int d = 0; d < v.Size(); d++)
3042  v(d) = p(d);
3043  }
3044  else
3045  {
3046  int d;
3047  for (d = 0; d < p.Size(); d++)
3048  v(d) = p(d);
3049  for ( ; d < v.Size(); d++)
3050  v(d) = 0.0;
3051  }
3052 }
3053 
3054 void Mesh::GetNodes(GridFunction &nodes) const
3055 {
3056  if (Nodes == NULL || Nodes->FESpace() != nodes.FESpace())
3057  {
3058  const int newSpaceDim = nodes.FESpace()->GetVDim();
3060  nodes.ProjectCoefficient(xyz);
3061  }
3062  else
3063  nodes = *Nodes;
3064 }
3065 
3067 {
3068  GridFunction *nodes = new GridFunction(nfes);
3069  SetNodalGridFunction(nodes, true);
3070 }
3071 
3072 void Mesh::SetNodalGridFunction(GridFunction *nodes, bool make_owner)
3073 {
3074  GetNodes(*nodes);
3075  NewNodes(*nodes, make_owner);
3076 }
3077 
3079 {
3080  return ((Nodes) ? Nodes->FESpace() : NULL);
3081 }
3082 
3084 {
3085  switch (Dim)
3086  {
3087  case 1: return GetNV();
3088  case 2: return GetNEdges();
3089  case 3: return GetNFaces();
3090  }
3091  return 0;
3092 }
3093 
3094 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
3095 static const char *fixed_or_not[] = { "fixed", "NOT FIXED" };
3096 #endif
3097 
3099 {
3100  int i, j, k, wo = 0, fo = 0, *vi = 0;
3101  double *v[4];
3102 
3103  if (Dim == 2 && spaceDim == 2)
3104  {
3105  DenseMatrix J(2, 2);
3106 
3107  for (i = 0; i < NumOfElements; i++)
3108  {
3109  if (Nodes == NULL)
3110  {
3111  vi = elements[i]->GetVertices();
3112  for (j = 0; j < 3; j++)
3113  v[j] = vertices[vi[j]]();
3114  for (j = 0; j < 2; j++)
3115  for (k = 0; k < 2; k++)
3116  J(j, k) = v[j+1][k] - v[0][k];
3117  }
3118  else
3119  {
3120  // only check the Jacobian at the center of the element
3121  GetElementJacobian(i, J);
3122  }
3123  if (J.Det() < 0.0)
3124  {
3125  if (fix_it)
3126  {
3127  switch (GetElementType(i))
3128  {
3129  case Element::TRIANGLE:
3130  mfem::Swap(vi[0], vi[1]);
3131  break;
3133  mfem::Swap(vi[1], vi[3]);
3134  break;
3135  }
3136  fo++;
3137  }
3138  wo++;
3139  }
3140  }
3141  }
3142 
3143  if (Dim == 3)
3144  {
3145  DenseMatrix J(3, 3);
3146 
3147  for (i = 0; i < NumOfElements; i++)
3148  {
3149  vi = elements[i]->GetVertices();
3150  switch (GetElementType(i))
3151  {
3152  case Element::TETRAHEDRON:
3153  if (Nodes == NULL)
3154  {
3155  for (j = 0; j < 4; j++)
3156  v[j] = vertices[vi[j]]();
3157  for (j = 0; j < 3; j++)
3158  for (k = 0; k < 3; k++)
3159  J(j, k) = v[j+1][k] - v[0][k];
3160  }
3161  else
3162  {
3163  // only check the Jacobian at the center of the element
3164  GetElementJacobian(i, J);
3165  }
3166  if (J.Det() < 0.0)
3167  {
3168  wo++;
3169  if (fix_it)
3170  {
3171  mfem::Swap(vi[0], vi[1]);
3172  fo++;
3173  }
3174  }
3175  break;
3176 
3177  case Element::HEXAHEDRON:
3178  // only check the Jacobian at the center of the element
3179  GetElementJacobian(i, J);
3180  if (J.Det() < 0.0)
3181  {
3182  wo++;
3183  if (fix_it)
3184  {
3185  // how?
3186  }
3187  }
3188  break;
3189  }
3190  }
3191  }
3192 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
3193  if (wo > 0)
3194  cout << "Elements with wrong orientation: " << wo << " / "
3195  << NumOfElements << " (" << fixed_or_not[(wo == fo) ? 0 : 1]
3196  << ")" << endl;
3197 #endif
3198 }
3199 
3200 int Mesh::GetTriOrientation(const int *base, const int *test)
3201 {
3202  int orient;
3203 
3204  if (test[0] == base[0])
3205  if (test[1] == base[1])
3206  orient = 0; // (0, 1, 2)
3207  else
3208  orient = 5; // (0, 2, 1)
3209  else if (test[0] == base[1])
3210  if (test[1] == base[0])
3211  orient = 1; // (1, 0, 2)
3212  else
3213  orient = 2; // (1, 2, 0)
3214  else // test[0] == base[2]
3215  if (test[1] == base[0])
3216  orient = 4; // (2, 0, 1)
3217  else
3218  orient = 3; // (2, 1, 0)
3219 
3220 #ifdef MFEM_DEBUG
3221  const int *aor = tri_orientations[orient];
3222  for (int j = 0; j < 3; j++)
3223  if (test[aor[j]] != base[j])
3224  mfem_error("Mesh::GetTriOrientation(...)");
3225 #endif
3226 
3227  return orient;
3228 }
3229 
3230 int Mesh::GetQuadOrientation(const int *base, const int *test)
3231 {
3232  int i;
3233 
3234  for (i = 0; i < 4; i++)
3235  if (test[i] == base[0])
3236  break;
3237 
3238 #ifdef MFEM_DEBUG
3239  int orient;
3240  if (test[(i+1)%4] == base[1])
3241  orient = 2*i;
3242  else
3243  orient = 2*i+1;
3244  const int *aor = quad_orientations[orient];
3245  for (int j = 0; j < 4; j++)
3246  if (test[aor[j]] != base[j])
3247  {
3248  cerr << "Mesh::GetQuadOrientation(...)" << endl;
3249  cerr << " base = [";
3250  for (int k = 0; k < 4; k++)
3251  cerr << " " << base[k];
3252  cerr << " ]\n test = [";
3253  for (int k = 0; k < 4; k++)
3254  cerr << " " << test[k];
3255  cerr << " ]" << endl;
3256  mfem_error();
3257  }
3258 #endif
3259 
3260  if (test[(i+1)%4] == base[1])
3261  return 2*i;
3262 
3263  return 2*i+1;
3264 }
3265 
3267 {
3268  int i, wo = 0;
3269 
3270  if (Dim == 2)
3271  {
3272  for (i = 0; i < NumOfBdrElements; i++)
3273  {
3274  if (faces_info[be_to_edge[i]].Elem2No < 0) // boundary face
3275  {
3276  int *bv = boundary[i]->GetVertices();
3277  int *fv = faces[be_to_edge[i]]->GetVertices();
3278  if (bv[0] != fv[0])
3279  {
3280  if (fix_it)
3281  mfem::Swap<int>(bv[0], bv[1]);
3282  wo++;
3283  }
3284  }
3285  }
3286  }
3287 
3288  if (Dim == 3)
3289  {
3290  int el, *bv, *ev;
3291  int v[4];
3292 
3293  for (i = 0; i < NumOfBdrElements; i++)
3294  {
3295  if (faces_info[be_to_face[i]].Elem2No < 0)
3296  { // boundary face
3297  bv = boundary[i]->GetVertices();
3298  el = faces_info[be_to_face[i]].Elem1No;
3299  ev = elements[el]->GetVertices();
3300  switch (GetElementType(el))
3301  {
3302  case Element::TETRAHEDRON:
3303  {
3304  int *fv = faces[be_to_face[i]]->GetVertices();
3305  int orientation; // orientation of the bdr. elem. w.r.t. the
3306  // corresponding face element (that's the base)
3307  orientation = GetTriOrientation(fv, bv);
3308  if (orientation % 2)
3309  {
3310  // wrong orientation -- swap vertices 0 and 1 so that
3311  // we don't change the marked edge: (0,1,2) -> (1,0,2)
3312  if (fix_it)
3313  mfem::Swap<int>(bv[0], bv[1]);
3314  wo++;
3315  }
3316  }
3317  break;
3318 
3319  case Element::HEXAHEDRON:
3320  {
3321  int lf = faces_info[be_to_face[i]].Elem1Inf/64;
3322  for (int j = 0; j < 4; j++)
3323  v[j] = ev[hex_faces[lf][j]];
3324  if (GetQuadOrientation(v, bv) % 2)
3325  {
3326  if (fix_it)
3327  mfem::Swap<int>(bv[0], bv[2]);
3328  wo++;
3329  }
3330  break;
3331  }
3332  }
3333  }
3334  }
3335  }
3336 // #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
3337 #ifdef MFEM_DEBUG
3338  if (wo > 0)
3339  cout << "Boundary elements with wrong orientation: " << wo << " / "
3340  << NumOfBdrElements << " (" << fixed_or_not[fix_it ? 0 : 1]
3341  << ")" << endl;
3342 #endif
3343 }
3344 
3346  const
3347 {
3348  if (el_to_edge)
3349  el_to_edge->GetRow(i, edges);
3350  else
3351  mfem_error("Mesh::GetElementEdges(...) element to edge table "
3352  "is not generated.");
3353 
3354  const int *v = elements[i]->GetVertices();
3355  const int ne = elements[i]->GetNEdges();
3356  cor.SetSize(ne);
3357  for (int j = 0; j < ne; j++)
3358  {
3359  const int *e = elements[i]->GetEdgeVertices(j);
3360  cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
3361  }
3362 }
3363 
3365  const
3366 {
3367  if (Dim == 2)
3368  {
3369  edges.SetSize(1);
3370  cor.SetSize(1);
3371  edges[0] = be_to_edge[i];
3372  const int *v = boundary[i]->GetVertices();
3373  cor[0] = (v[0] < v[1]) ? (1) : (-1);
3374  }
3375  else if (Dim == 3)
3376  {
3377  if (bel_to_edge)
3378  bel_to_edge->GetRow(i, edges);
3379  else
3380  mfem_error("Mesh::GetBdrElementEdges(...)");
3381 
3382  const int *v = boundary[i]->GetVertices();
3383  const int ne = boundary[i]->GetNEdges();
3384  cor.SetSize(ne);
3385  for (int j = 0; j < ne; j++)
3386  {
3387  const int *e = boundary[i]->GetEdgeVertices(j);
3388  cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
3389  }
3390  }
3391 }
3392 
3393 void Mesh::GetFaceEdges(int i, Array<int> &edges, Array<int> &o) const
3394 {
3395  if (Dim == 2)
3396  {
3397  edges.SetSize(1);
3398  edges[0] = i;
3399  o.SetSize(1);
3400  const int *v = faces[i]->GetVertices();
3401  o[0] = (v[0] < v[1]) ? (1) : (-1);
3402  }
3403 
3404  if (Dim != 3)
3405  return;
3406 
3407  GetFaceEdgeTable(); // generate face_edge Table (if not generated)
3408 
3409  face_edge->GetRow(i, edges);
3410 
3411  const int *v = faces[i]->GetVertices();
3412  const int ne = faces[i]->GetNEdges();
3413  o.SetSize(ne);
3414  for (int j = 0; j < ne; j++)
3415  {
3416  const int *e = faces[i]->GetEdgeVertices(j);
3417  o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
3418  }
3419 }
3420 
3421 void Mesh::GetEdgeVertices(int i, Array<int> &vert) const
3422 {
3423  // the two vertices are sorted: vert[0] < vert[1]
3424  // this is consistent with the global edge orientation
3425  GetEdgeVertexTable(); // generate edge_vertex Table (if not generated)
3426  edge_vertex->GetRow(i, vert);
3427 }
3428 
3430 {
3431  if (face_edge)
3432  return face_edge;
3433 
3434  if (Dim != 3)
3435  return NULL;
3436 
3437 #ifdef MFEM_DEBUG
3438  if (faces.Size() != NumOfFaces)
3439  mfem_error("Mesh::GetFaceEdgeTable : faces were not generated!");
3440 #endif
3441 
3442  DSTable v_to_v(NumOfVertices);
3443  GetVertexToVertexTable(v_to_v);
3444 
3445  face_edge = new Table;
3446  GetElementArrayEdgeTable(faces, v_to_v, *face_edge);
3447 
3448  return(face_edge);
3449 }
3450 
3452 {
3453  if (edge_vertex)
3454  return edge_vertex;
3455 
3456  DSTable v_to_v(NumOfVertices);
3457  GetVertexToVertexTable(v_to_v);
3458 
3459  int nedges = v_to_v.NumberOfEntries();
3460  edge_vertex = new Table(nedges, 2);
3461  for (int i = 0; i < NumOfVertices; i++)
3462  {
3463  for (DSTable::RowIterator it(v_to_v, i); !it; ++it)
3464  {
3465  int j = it.Index();
3466  edge_vertex->Push(j, i);
3467  edge_vertex->Push(j, it.Column());
3468  }
3469  }
3470  edge_vertex->Finalize();
3471 
3472  return edge_vertex;
3473 }
3474 
3476 {
3477  int i, j, nv, *v;
3478 
3479  Table *vert_elem = new Table;
3480 
3481  vert_elem->MakeI(NumOfVertices);
3482 
3483  for (i = 0; i < NumOfElements; i++)
3484  {
3485  nv = elements[i]->GetNVertices();
3486  v = elements[i]->GetVertices();
3487  for (j = 0; j < nv; j++)
3488  vert_elem->AddAColumnInRow(v[j]);
3489  }
3490 
3491  vert_elem->MakeJ();
3492 
3493  for (i = 0; i < NumOfElements; i++)
3494  {
3495  nv = elements[i]->GetNVertices();
3496  v = elements[i]->GetVertices();
3497  for (j = 0; j < nv; j++)
3498  vert_elem->AddConnection(v[j], i);
3499  }
3500 
3501  vert_elem->ShiftUpI();
3502 
3503  return vert_elem;
3504 }
3505 
3507 {
3508  Table *face_elem = new Table;
3509 
3510  face_elem->MakeI(faces_info.Size());
3511 
3512  for (int i = 0; i < faces_info.Size(); i++)
3513  {
3514  if (faces_info[i].Elem2No >= 0)
3515  face_elem->AddColumnsInRow(i, 2);
3516  else
3517  face_elem->AddAColumnInRow(i);
3518  }
3519 
3520  face_elem->MakeJ();
3521 
3522  for (int i = 0; i < faces_info.Size(); i++)
3523  {
3524  face_elem->AddConnection(i, faces_info[i].Elem1No);
3525  if (faces_info[i].Elem2No >= 0)
3526  face_elem->AddConnection(i, faces_info[i].Elem2No);
3527  }
3528 
3529  face_elem->ShiftUpI();
3530 
3531  return face_elem;
3532 }
3533 
3535  const
3536 {
3537  int n, j;
3538 
3539  if (el_to_face)
3540  el_to_face->GetRow(i, fcs);
3541  else
3542  mfem_error("Mesh::GetElementFaces(...) : el_to_face not generated.");
3543 
3544  n = fcs.Size();
3545  cor.SetSize(n);
3546  for (j = 0; j < n; j++)
3547  if (faces_info[fcs[j]].Elem1No == i)
3548  cor[j] = faces_info[fcs[j]].Elem1Inf % 64;
3549 #ifdef MFEM_DEBUG
3550  else if (faces_info[fcs[j]].Elem2No == i)
3551  cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
3552  else
3553  mfem_error("Mesh::GetElementFaces(...) : 2");
3554 #else
3555  else
3556  cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
3557 #endif
3558 }
3559 
3560 void Mesh::GetBdrElementFace(int i, int *f, int *o) const
3561 {
3562  const int *bv, *fv;
3563 
3564  if (State == Mesh::TWO_LEVEL_COARSE)
3565  {
3566  // the coarse level 'be_to_face' and 'faces' are destroyed
3567  mfem_error("Mesh::GetBdrElementFace (...)");
3568  }
3569 
3570  *f = be_to_face[i];
3571  bv = boundary[i]->GetVertices();
3572  fv = faces[be_to_face[i]]->GetVertices();
3573 
3574  // find the orientation of the bdr. elem. w.r.t.
3575  // the corresponding face element (that's the base)
3576  switch (GetBdrElementType(i))
3577  {
3578  case Element::TRIANGLE:
3579  *o = GetTriOrientation(fv, bv);
3580  break;
3582  *o = GetQuadOrientation(fv, bv);
3583  break;
3584  default:
3585  mfem_error("Mesh::GetBdrElementFace(...) 2");
3586  }
3587 }
3588 
3590 {
3591  // Here, we assume all faces are of the same type
3592  switch (GetElementType(0))
3593  {
3594  case Element::SEGMENT:
3595  return Geometry::POINT;
3596 
3597  case Element::TRIANGLE:
3599  return Geometry::SEGMENT; // in 2D 'face' is an edge
3600 
3601  case Element::TETRAHEDRON:
3602  return Geometry::TRIANGLE;
3603  case Element::HEXAHEDRON:
3604  return Geometry::SQUARE;
3605  default:
3606  mfem_error("Mesh::GetFaceBaseGeometry(...) #1");
3607  }
3608  return(-1);
3609 #if 0
3610  if (faces[i] == NULL)
3611  switch (GetElementType(faces_info[i].Elem1No))
3612  {
3613  case Element::TETRAHEDRON:
3614  return Geometry::TRIANGLE;
3615  case Element::HEXAHEDRON:
3616  return Geometry::SQUARE;
3617  default:
3618  mfem_error("Mesh::GetFaceBaseGeometry(...) #2");
3619  }
3620  else
3621  return faces[i]->GetGeometryType();
3622 #endif
3623 }
3624 
3626 {
3627  if (Dim == 2)
3628  return be_to_edge[i];
3629  return be_to_face[i];
3630 }
3631 
3632 int Mesh::GetElementType(int i) const
3633 {
3634  Element *El = elements[i];
3635  int t = El->GetType();
3636 
3637  while (1)
3638  if (t == Element::BISECTED ||
3639  t == Element::QUADRISECTED ||
3640  t == Element::OCTASECTED)
3641  t = (El = ((RefinedElement *) El)->IAm())->GetType();
3642  else
3643  break;
3644  return t;
3645 }
3646 
3647 int Mesh::GetBdrElementType(int i) const
3648 {
3649  Element *El = boundary[i];
3650  int t = El->GetType();
3651 
3652  while (1)
3653  if (t == Element::BISECTED || t == Element::QUADRISECTED)
3654  t = (El = ((RefinedElement *) El)->IAm())->GetType();
3655  else
3656  break;
3657  return t;
3658 }
3659 
3660 void Mesh::GetPointMatrix(int i, DenseMatrix &pointmat) const
3661 {
3662  int k, j, nv;
3663  const int *v;
3664 
3665  v = elements[i]->GetVertices();
3666  nv = elements[i]->GetNVertices();
3667 
3668  pointmat.SetSize(spaceDim, nv);
3669  for (k = 0; k < spaceDim; k++)
3670  for (j = 0; j < nv; j++)
3671  pointmat(k, j) = vertices[v[j]](k);
3672 }
3673 
3674 void Mesh::GetBdrPointMatrix(int i,DenseMatrix &pointmat) const
3675 {
3676  int k, j, nv;
3677  const int *v;
3678 
3679  v = boundary[i]->GetVertices();
3680  nv = boundary[i]->GetNVertices();
3681 
3682  pointmat.SetSize(spaceDim, nv);
3683  for (k = 0; k < spaceDim; k++)
3684  for (j = 0; j < nv; j++)
3685  pointmat(k, j) = vertices[v[j]](k);
3686 }
3687 
3688 double Mesh::GetLength(int i, int j) const
3689 {
3690  const double *vi = vertices[i]();
3691  const double *vj = vertices[j]();
3692  double length = 0.;
3693 
3694  for (int k = 0; k < spaceDim; k++)
3695  length += (vi[k]-vj[k])*(vi[k]-vj[k]);
3696 
3697  return sqrt(length);
3698 }
3699 
3700 // static method
3702  const DSTable &v_to_v, Table &el_to_edge)
3703 {
3704  el_to_edge.MakeI(elem_array.Size());
3705  for (int i = 0; i < elem_array.Size(); i++)
3706  {
3707  el_to_edge.AddColumnsInRow(i, elem_array[i]->GetNEdges());
3708  }
3709  el_to_edge.MakeJ();
3710  for (int i = 0; i < elem_array.Size(); i++)
3711  {
3712  const int *v = elem_array[i]->GetVertices();
3713  const int ne = elem_array[i]->GetNEdges();
3714  for (int j = 0; j < ne; j++)
3715  {
3716  const int *e = elem_array[i]->GetEdgeVertices(j);
3717  el_to_edge.AddConnection(i, v_to_v(v[e[0]], v[e[1]]));
3718  }
3719  }
3720  el_to_edge.ShiftUpI();
3721 }
3722 
3724 {
3725  if (edge_vertex)
3726  {
3727  for (int i = 0; i < edge_vertex->Size(); i++)
3728  {
3729  const int *v = edge_vertex->GetRow(i);
3730  v_to_v.Push(v[0], v[1]);
3731  }
3732  }
3733  else
3734  {
3735  for (int i = 0; i < NumOfElements; i++)
3736  {
3737  const int *v = elements[i]->GetVertices();
3738  const int ne = elements[i]->GetNEdges();
3739  for (int j = 0; j < ne; j++)
3740  {
3741  const int *e = elements[i]->GetEdgeVertices(j);
3742  v_to_v.Push(v[e[0]], v[e[1]]);
3743  }
3744  }
3745  }
3746 }
3747 
3749 {
3750  int i, NumberOfEdges;
3751 
3752  DSTable v_to_v(NumOfVertices);
3753  GetVertexToVertexTable(v_to_v);
3754 
3755  NumberOfEdges = v_to_v.NumberOfEntries();
3756 
3757  // Fill the element to edge table
3758  GetElementArrayEdgeTable(elements, v_to_v, e_to_f);
3759 
3760  if (Dim == 2)
3761  {
3762  // Initialize the indeces for the boundary elements.
3763  be_to_f.SetSize(NumOfBdrElements);
3764  for (i = 0; i < NumOfBdrElements; i++)
3765  {
3766  const int *v = boundary[i]->GetVertices();
3767  be_to_f[i] = v_to_v(v[0], v[1]);
3768  }
3769  }
3770  else if (Dim == 3)
3771  {
3772  if (bel_to_edge == NULL)
3773  bel_to_edge = new Table;
3774  GetElementArrayEdgeTable(boundary, v_to_v, *bel_to_edge);
3775  }
3776  else
3777  mfem_error("1D GetElementToEdgeTable is not yet implemented.");
3778 
3779  // Return the number of edges
3780  return NumberOfEdges;
3781 }
3782 
3784 {
3785  if (el_to_el)
3786  return *el_to_el;
3787 
3788  if (Dim < 3)
3789  {
3790  Table *f_el;
3791  if (Dim == 1)
3792  {
3793  f_el = GetVertexToElementTable();
3794  }
3795  else
3796  {
3797  f_el = new Table;
3798  Transpose(ElementToEdgeTable(), *f_el);
3799  }
3800 
3801  el_to_el = new Table;
3802  el_to_el->MakeI(NumOfElements);
3803  for (int i = 0; i < f_el->Size(); i++)
3804  {
3805  if (f_el->RowSize(i) > 1)
3806  {
3807  const int *el = f_el->GetRow(i);
3808  el_to_el->AddAColumnInRow(el[0]);
3809  el_to_el->AddAColumnInRow(el[1]);
3810  }
3811  }
3812  el_to_el->MakeJ();
3813  for (int i = 0; i < f_el->Size(); i++)
3814  {
3815  if (f_el->RowSize(i) > 1)
3816  {
3817  const int *el = f_el->GetRow(i);
3818  el_to_el->AddConnection(el[0], el[1]);
3819  el_to_el->AddConnection(el[1], el[0]);
3820  }
3821  }
3822  el_to_el->ShiftUpI();
3823 
3824  delete f_el;
3825  }
3826  else
3827  {
3828  el_to_el = new Table(NumOfElements, 6); // 6 is the max. # of faces
3829 
3830 #ifdef MFEM_DEBUG
3831  if (faces_info.Size() != NumOfFaces)
3832  mfem_error("Mesh::ElementToElementTable : faces were not generated!");
3833 #endif
3834 
3835  for (int i = 0; i < faces_info.Size(); i++)
3836  if (faces_info[i].Elem2No >= 0)
3837  {
3838  el_to_el->Push(faces_info[i].Elem1No, faces_info[i].Elem2No);
3839  el_to_el->Push(faces_info[i].Elem2No, faces_info[i].Elem1No);
3840  }
3841 
3842  el_to_el->Finalize();
3843  }
3844 
3845  return *el_to_el;
3846 }
3847 
3849 {
3850  if (el_to_face == NULL)
3851  mfem_error("Mesh::ElementToFaceTable()");
3852  return *el_to_face;
3853 }
3854 
3856 {
3857  if (el_to_edge == NULL)
3858  mfem_error("Mesh::ElementToEdgeTable()");
3859  return *el_to_edge;
3860 }
3861 
3862 void Mesh::AddPointFaceElement(int lf, int gf, int el)
3863 {
3864  if (faces_info[gf].Elem1No == -1) // this will be elem1
3865  {
3866  // faces[gf] = new Point(&gf);
3867  faces_info[gf].Elem1No = el;
3868  faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
3869  faces_info[gf].Elem2No = -1; // in case there's no other side
3870  faces_info[gf].Elem2Inf = -1; // face is not shared
3871  }
3872  else // this will be elem2
3873  {
3874 #ifdef MFEM_DEBUG
3875  // int *v = faces[gf]->GetVertices();
3876  // if (v[0] != gf)
3877  // mfem_error("Mesh::AddPointFaceElement(...)");
3878 #endif
3879  faces_info[gf].Elem2No = el;
3880  faces_info[gf].Elem2Inf = 64 * lf + 1;
3881  }
3882 }
3883 
3884 void Mesh::AddSegmentFaceElement(int lf, int gf, int el, int v0, int v1)
3885 {
3886  if (faces[gf] == NULL) // this will be elem1
3887  {
3888  faces[gf] = new Segment(v0, v1);
3889  faces_info[gf].Elem1No = el;
3890  faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
3891  faces_info[gf].Elem2No = -1; // in case there's no other side
3892  faces_info[gf].Elem2Inf = -1; // face is not shared
3893  }
3894  else // this will be elem2
3895  {
3896 #ifdef MFEM_DEBUG
3897  int *v = faces[gf]->GetVertices();
3898  if (v[1] != v0 || v[0] != v1)
3899  mfem_error("Mesh::AddSegmentFaceElement(...)");
3900 #endif
3901  faces_info[gf].Elem2No = el;
3902  faces_info[gf].Elem2Inf = 64 * lf + 1;
3903  }
3904 }
3905 
3906 void Mesh::AddTriangleFaceElement(int lf, int gf, int el,
3907  int v0, int v1, int v2)
3908 {
3909  if (faces[gf] == NULL) // this will be elem1
3910  {
3911  faces[gf] = new Triangle(v0, v1, v2);
3912  faces_info[gf].Elem1No = el;
3913  faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
3914  faces_info[gf].Elem2No = -1; // in case there's no other side
3915  faces_info[gf].Elem2Inf = -1; // face is not shared
3916  }
3917  else // this will be elem2
3918  {
3919  int orientation, vv[3] = { v0, v1, v2 };
3920  orientation = GetTriOrientation(faces[gf]->GetVertices(), vv);
3921 #ifdef MFEM_DEBUG
3922  if (orientation % 2 == 0)
3923  mfem_error("Mesh::AddTriangleFaceElement(...)");
3924 #endif
3925  faces_info[gf].Elem2No = el;
3926  faces_info[gf].Elem2Inf = 64 * lf + orientation;
3927  }
3928 }
3929 
3930 void Mesh::AddQuadFaceElement(int lf, int gf, int el,
3931  int v0, int v1, int v2, int v3)
3932 {
3933  if (faces_info[gf].Elem1No < 0) // this will be elem1
3934  {
3935  faces[gf] = new Quadrilateral(v0, v1, v2, v3);
3936  faces_info[gf].Elem1No = el;
3937  faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
3938  faces_info[gf].Elem2No = -1; // in case there's no other side
3939  faces_info[gf].Elem2Inf = -1; // face is not shared
3940  }
3941  else // this will be elem2
3942  {
3943  int vv[4] = { v0, v1, v2, v3 };
3944  int oo = GetQuadOrientation(faces[gf]->GetVertices(), vv);
3945 #ifdef MFEM_DEBUG
3946  if (oo % 2 == 0)
3947  mfem_error("Mesh::AddQuadFaceElement(...)");
3948 #endif
3949  faces_info[gf].Elem2No = el;
3950  faces_info[gf].Elem2Inf = 64 * lf + oo;
3951  }
3952 }
3953 
3955 {
3956  int i, nfaces;
3957 
3958  nfaces = (Dim == 1) ? NumOfVertices : ((Dim == 2) ? NumOfEdges : NumOfFaces);
3959 
3960  for (i = 0; i < faces.Size(); i++)
3961  FreeElement(faces[i]);
3962 
3963  // (re)generate the interior faces and the info for them
3964  faces.SetSize(nfaces);
3965  faces_info.SetSize(nfaces);
3966  for (i = 0; i < nfaces; i++)
3967  {
3968  faces[i] = NULL;
3969  faces_info[i].Elem1No = -1;
3970  }
3971  for (i = 0; i < NumOfElements; i++)
3972  {
3973  const int *v = elements[i]->GetVertices();
3974  const int *ef;
3975  if (Dim == 1)
3976  {
3977  AddPointFaceElement(0, v[0], i);
3978  AddPointFaceElement(1, v[1], i);
3979  }
3980  else if (Dim == 2)
3981  {
3982  ef = el_to_edge->GetRow(i);
3983  const int ne = elements[i]->GetNEdges();
3984  for (int j = 0; j < ne; j++)
3985  {
3986  const int *e = elements[i]->GetEdgeVertices(j);
3987  AddSegmentFaceElement(j, ef[j], i, v[e[0]], v[e[1]]);
3988  }
3989  }
3990  else
3991  {
3992  ef = el_to_face->GetRow(i);
3993  switch (GetElementType(i))
3994  {
3995  case Element::TETRAHEDRON:
3996  for (int j = 0; j < 4; j++)
3997  {
3998  const int *fv = tet_faces[j];
3999  AddTriangleFaceElement(j, ef[j], i,
4000  v[fv[0]], v[fv[1]], v[fv[2]]);
4001  }
4002  break;
4003  case Element::HEXAHEDRON:
4004  for (int j = 0; j < 6; j++)
4005  {
4006  const int *fv = hex_faces[j];
4007  AddQuadFaceElement(j, ef[j], i,
4008  v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
4009  }
4010  break;
4011 #ifdef MFEM_DEBUG
4012  default:
4013  MFEM_ABORT("Unexpected type of Element.");
4014 #endif
4015  }
4016  }
4017  }
4018 }
4019 
4021 {
4022  STable3D *faces_tbl = new STable3D(NumOfVertices);
4023  for (int i = 0; i < NumOfElements; i++)
4024  {
4025  const int *v = elements[i]->GetVertices();
4026  switch (GetElementType(i))
4027  {
4028  case Element::TETRAHEDRON:
4029  for (int j = 0; j < 4; j++)
4030  {
4031  const int *fv = tet_faces[j];
4032  faces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
4033  }
4034  break;
4035  case Element::HEXAHEDRON:
4036  // find the face by the vertices with the smallest 3 numbers
4037  // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
4038  for (int j = 0; j < 6; j++)
4039  {
4040  const int *fv = hex_faces[j];
4041  faces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
4042  }
4043  break;
4044 #ifdef MFEM_DEBUG
4045  default:
4046  MFEM_ABORT("Unexpected type of Element.");
4047 #endif
4048  }
4049  }
4050  return faces_tbl;
4051 }
4052 
4054 {
4055  int i, *v;
4056  STable3D *faces_tbl;
4057 
4058  if (el_to_face != NULL)
4059  delete el_to_face;
4060  el_to_face = new Table(NumOfElements, 6); // must be 6 for hexahedra
4061  faces_tbl = new STable3D(NumOfVertices);
4062  for (i = 0; i < NumOfElements; i++)
4063  {
4064  v = elements[i]->GetVertices();
4065  switch (GetElementType(i))
4066  {
4067  case Element::TETRAHEDRON:
4068  for (int j = 0; j < 4; j++)
4069  {
4070  const int *fv = tet_faces[j];
4071  el_to_face->Push(
4072  i, faces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]));
4073  }
4074  break;
4075  case Element::HEXAHEDRON:
4076  // find the face by the vertices with the smallest 3 numbers
4077  // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
4078  for (int j = 0; j < 6; j++)
4079  {
4080  const int *fv = hex_faces[j];
4081  el_to_face->Push(
4082  i, faces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
4083  }
4084  break;
4085 #ifdef MFEM_DEBUG
4086  default:
4087  MFEM_ABORT("Unexpected type of Element.");
4088 #endif
4089  }
4090  }
4091  el_to_face->Finalize();
4092  NumOfFaces = faces_tbl->NumberOfElements();
4093  be_to_face.SetSize(NumOfBdrElements);
4094  for (i = 0; i < NumOfBdrElements; i++)
4095  {
4096  v = boundary[i]->GetVertices();
4097  switch (GetBdrElementType(i))
4098  {
4099  case Element::TRIANGLE:
4100  be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
4101  break;
4103  be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
4104  break;
4105 #ifdef MFEM_DEBUG
4106  default:
4107  MFEM_ABORT("Unexpected type of boundary Element.");
4108 #endif
4109  }
4110  }
4111 
4112  if (ret_ftbl)
4113  return faces_tbl;
4114  delete faces_tbl;
4115  return NULL;
4116 }
4117 
4119 {
4120  int *v;
4121 
4122  if (Dim != 3 || !(meshgen & 1))
4123  return;
4124 
4125  DSTable *old_v_to_v = NULL;
4126  Table *old_elem_vert = NULL;
4127 
4128  if (Nodes)
4129  PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
4130 
4131  DeleteCoarseTables();
4132 
4133  for (int i = 0; i < NumOfElements; i++)
4134  if (GetElementType(i) == Element::TETRAHEDRON)
4135  {
4136  v = elements[i]->GetVertices();
4137 
4138  Rotate3(v[0], v[1], v[2]);
4139  if (v[0] < v[3])
4140  Rotate3(v[1], v[2], v[3]);
4141  else
4142  ShiftL2R(v[0], v[1], v[3]);
4143  }
4144 
4145  for (int i = 0; i < NumOfBdrElements; i++)
4146  if (GetBdrElementType(i) == Element::TRIANGLE)
4147  {
4148  v = boundary[i]->GetVertices();
4149 
4150  Rotate3(v[0], v[1], v[2]);
4151  }
4152 
4153  if (!Nodes)
4154  {
4155  GetElementToFaceTable();
4156  GenerateFaces();
4157  if (el_to_edge)
4158  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
4159  }
4160  else
4161  {
4162  DoNodeReorder(old_v_to_v, old_elem_vert);
4163  delete old_elem_vert;
4164  delete old_v_to_v;
4165  }
4166 }
4167 
4168 #ifdef MFEM_USE_MPI
4169 // auxiliary function for qsort
4170 static int mfem_less(const void *x, const void *y)
4171 {
4172  if (*(int*)x < *(int*)y)
4173  return 1;
4174  if (*(int*)x > *(int*)y)
4175  return -1;
4176  return 0;
4177 }
4178 #ifndef MFEM_USE_METIS_5
4179 // METIS 4 prototypes
4180 typedef int idxtype;
4181 extern "C" {
4183  int*, int*, int*, int*, int*, idxtype*);
4185  int*, int*, int*, int*, int*, idxtype*);
4187  int*, int*, int*, int*, int*, idxtype*);
4188 }
4189 #else
4190 // METIS 5 prototypes
4191 typedef int idx_t;
4192 typedef float real_t;
4193 extern "C" {
4195  idx_t *nvtxs, idx_t *ncon, idx_t *xadj,
4196  idx_t *adjncy, idx_t *vwgt, idx_t *vsize, idx_t *adjwgt,
4197  idx_t *nparts, real_t *tpwgts, real_t *ubvec, idx_t *options,
4198  idx_t *edgecut, idx_t *part);
4199  int METIS_PartGraphKway(
4200  idx_t *nvtxs, idx_t *ncon, idx_t *xadj,
4201  idx_t *adjncy, idx_t *vwgt, idx_t *vsize, idx_t *adjwgt,
4202  idx_t *nparts, real_t *tpwgts, real_t *ubvec, idx_t *options,
4203  idx_t *edgecut, idx_t *part);
4204  int METIS_SetDefaultOptions(idx_t *options);
4205 }
4206 #endif
4207 #endif
4208 
4210 {
4211  int *partitioning;
4212  double pmin[3] = { numeric_limits<double>::infinity(),
4213  numeric_limits<double>::infinity(),
4214  numeric_limits<double>::infinity()};
4215  double pmax[3] = { -numeric_limits<double>::infinity(),
4216  -numeric_limits<double>::infinity(),
4217  -numeric_limits<double>::infinity()};
4218  // find a bounding box using the vertices
4219  for (int vi = 0; vi < NumOfVertices; vi++)
4220  {
4221  const double *p = vertices[vi]();
4222  for (int i = 0; i < spaceDim; i++)
4223  {
4224  if (p[i] < pmin[i]) pmin[i] = p[i];
4225  if (p[i] > pmax[i]) pmax[i] = p[i];
4226  }
4227  }
4228 
4229  partitioning = new int[NumOfElements];
4230 
4231  // determine the partitioning using the centers of the elements
4232  double ppt[3];
4233  Vector pt(ppt, spaceDim);
4234  for (int el = 0; el < NumOfElements; el++)
4235  {
4236  GetElementTransformation(el)->Transform(
4237  Geometries.GetCenter(GetElementBaseGeometry(el)), pt);
4238  int part = 0;
4239  for (int i = spaceDim-1; i >= 0; i--)
4240  {
4241  int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
4242  if (idx < 0) idx = 0;
4243  if (idx >= nxyz[i]) idx = nxyz[i]-1;
4244  part = part * nxyz[i] + idx;
4245  }
4246  partitioning[el] = part;
4247  }
4248 
4249  return partitioning;
4250 }
4251 
4252 int *Mesh::GeneratePartitioning(int nparts, int part_method)
4253 {
4254 #ifdef MFEM_USE_MPI
4255  int i, *partitioning;
4256 
4257  ElementToElementTable();
4258 
4259  partitioning = new int[NumOfElements];
4260 
4261  if (nparts == 1)
4262  {
4263  for (i = 0; i < NumOfElements; i++)
4264  partitioning[i] = 0;
4265  }
4266  else
4267  {
4268  int *I, *J, n;
4269 #ifndef MFEM_USE_METIS_5
4270  int wgtflag = 0;
4271  int numflag = 0;
4272  int options[5];
4273 #else
4274  int ncon = 1;
4275  int err;
4276  int options[40];
4277 #endif
4278  int edgecut;
4279 
4280  n = NumOfElements;
4281  I = el_to_el->GetI();
4282  J = el_to_el->GetJ();
4283 #ifndef MFEM_USE_METIS_5
4284  options[0] = 0;
4285 #else
4286  METIS_SetDefaultOptions(options);
4287  options[10] = 1; // set METIS_OPTION_CONTIG
4288 #endif
4289 
4290  // Sort the neighbor lists
4291  if (part_method >= 0 && part_method <= 2)
4292  for (i = 0; i < n; i++)
4293  qsort(&J[I[i]], I[i+1]-I[i], sizeof(int), &mfem_less);
4294 
4295  // This function should be used to partition a graph into a small
4296  // number of partitions (less than 8).
4297  if (part_method == 0 || part_method == 3)
4298  {
4299 #ifndef MFEM_USE_METIS_5
4301  (idxtype *) I,
4302  (idxtype *) J,
4303  (idxtype *) NULL,
4304  (idxtype *) NULL,
4305  &wgtflag,
4306  &numflag,
4307  &nparts,
4308  options,
4309  &edgecut,
4310  (idxtype *) partitioning);
4311 #else
4312  err = METIS_PartGraphRecursive(&n,
4313  &ncon,
4314  I,
4315  J,
4316  (idx_t *) NULL,
4317  (idx_t *) NULL,
4318  (idx_t *) NULL,
4319  &nparts,
4320  (real_t *) NULL,
4321  (real_t *) NULL,
4322  options,
4323  &edgecut,
4324  partitioning);
4325  if (err != 1)
4326  mfem_error("Mesh::GeneratePartitioning: "
4327  " error in METIS_PartGraphRecursive!");
4328 #endif
4329  }
4330 
4331  // This function should be used to partition a graph into a large
4332  // number of partitions (greater than 8).
4333  if (part_method == 1 || part_method == 4)
4334  {
4335 #ifndef MFEM_USE_METIS_5
4337  (idxtype *) I,
4338  (idxtype *) J,
4339  (idxtype *) NULL,
4340  (idxtype *) NULL,
4341  &wgtflag,
4342  &numflag,
4343  &nparts,
4344  options,
4345  &edgecut,
4346  (idxtype *) partitioning);
4347 #else
4348  err = METIS_PartGraphKway(&n,
4349  &ncon,
4350  I,
4351  J,
4352  (idx_t *) NULL,
4353  (idx_t *) NULL,
4354  (idx_t *) NULL,
4355  &nparts,
4356  (real_t *) NULL,
4357  (real_t *) NULL,
4358  options,
4359  &edgecut,
4360  partitioning);
4361  if (err != 1)
4362  mfem_error("Mesh::GeneratePartitioning: "
4363  " error in METIS_PartGraphKway!");
4364 #endif
4365  }
4366 
4367  // The objective of this partitioning is to minimize the total
4368  // communication volume
4369  if (part_method == 2 || part_method == 5)
4370  {
4371 #ifndef MFEM_USE_METIS_5
4373  (idxtype *) I,
4374  (idxtype *) J,
4375  (idxtype *) NULL,
4376  (idxtype *) NULL,
4377  &wgtflag,
4378  &numflag,
4379  &nparts,
4380  options,
4381  &edgecut,
4382  (idxtype *) partitioning);
4383 #else
4384  options[1] = 1; // set METIS_OPTION_OBJTYPE to METIS_OBJTYPE_VOL
4385  err = METIS_PartGraphKway(&n,
4386  &ncon,
4387  I,
4388  J,
4389  (idx_t *) NULL,
4390  (idx_t *) NULL,
4391  (idx_t *) NULL,
4392  &nparts,
4393  (real_t *) NULL,
4394  (real_t *) NULL,
4395  options,
4396  &edgecut,
4397  partitioning);
4398  if (err != 1)
4399  mfem_error("Mesh::GeneratePartitioning: "
4400  " error in METIS_PartGraphKway!");
4401 #endif
4402  }
4403 
4404 #ifdef MFEM_DEBUG
4405  cout << "Mesh::GeneratePartitioning(...): edgecut = "
4406  << edgecut << endl;
4407 #endif
4408  }
4409 
4410  if (el_to_el)
4411  delete el_to_el;
4412  el_to_el = NULL;
4413 
4414  // Check for empty partitionings (a "feature" in METIS)
4415  {
4416  Array< Pair<int,int> > psize(nparts);
4417  for (i = 0; i < nparts; i++)
4418  {
4419  psize[i].one = 0;
4420  psize[i].two = i;
4421  }
4422 
4423  for (i = 0; i < NumOfElements; i++)
4424  psize[partitioning[i]].one++;
4425 
4426  int empty_parts = 0;
4427  for (i = 0; i < nparts; i++)
4428  if (psize[i].one == 0)
4429  empty_parts++;
4430 
4431  // This code just split the largest partitionings in two.
4432  // Do we need to replace it with something better?
4433  if (empty_parts)
4434  {
4435  cerr << "Mesh::GeneratePartitioning returned " << empty_parts
4436  << " empty parts!" << endl;
4437 
4438  SortPairs<int,int>(psize, nparts);
4439 
4440  for (i = nparts-1; i > nparts-1-empty_parts; i--)
4441  psize[i].one /= 2;
4442 
4443  for (int j = 0; j < NumOfElements; j++)
4444  for (i = nparts-1; i > nparts-1-empty_parts; i--)
4445  if (psize[i].one == 0 || partitioning[j] != psize[i].two)
4446  continue;
4447  else
4448  {
4449  partitioning[j] = psize[nparts-1-i].two;
4450  psize[i].one--;
4451  }
4452  }
4453  }
4454 
4455  return partitioning;
4456 
4457 #else
4458 
4459  mfem_error("Mesh::GeneratePartitioning(...): "
4460  "MFEM was compiled without Metis.");
4461 
4462  return NULL;
4463 
4464 #endif
4465 }
4466 
4467 /* required: 0 <= partitioning[i] < num_part */
4469  const Array<int> &partitioning,
4470  Array<int> &component,
4471  Array<int> &num_comp)
4472 {
4473  int i, j, k;
4474  int num_elem, *i_elem_elem, *j_elem_elem;
4475 
4476  num_elem = elem_elem.Size();
4477  i_elem_elem = elem_elem.GetI();
4478  j_elem_elem = elem_elem.GetJ();
4479 
4480  component.SetSize(num_elem);
4481 
4482  Array<int> elem_stack(num_elem);
4483  int stack_p, stack_top_p, elem;
4484  int num_part;
4485 
4486  num_part = -1;
4487  for (i = 0; i < num_elem; i++)
4488  {
4489  if (partitioning[i] > num_part)
4490  num_part = partitioning[i];
4491  component[i] = -1;
4492  }
4493  num_part++;
4494 
4495  num_comp.SetSize(num_part);
4496  for (i = 0; i < num_part; i++)
4497  num_comp[i] = 0;
4498 
4499  stack_p = 0;
4500  stack_top_p = 0; // points to the first unused element in the stack
4501  for (elem = 0; elem < num_elem; elem++)
4502  {
4503  if (component[elem] >= 0)
4504  continue;
4505 
4506  component[elem] = num_comp[partitioning[elem]]++;
4507 
4508  elem_stack[stack_top_p++] = elem;
4509 
4510  for ( ; stack_p < stack_top_p; stack_p++)
4511  {
4512  i = elem_stack[stack_p];
4513  for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
4514  {
4515  k = j_elem_elem[j];
4516  if (partitioning[k] == partitioning[i])
4517  {
4518  if (component[k] < 0)
4519  {
4520  component[k] = component[i];
4521  elem_stack[stack_top_p++] = k;
4522  }
4523  else if (component[k] != component[i])
4524  {
4525  mfem_error("FindPartitioningComponents");
4526  }
4527  }
4528  }
4529  }
4530  }
4531 }
4532 
4533 void Mesh::CheckPartitioning(int *partitioning)
4534 {
4535  int i, n_empty, n_mcomp;
4536  Array<int> component, num_comp;
4537  const Array<int> _partitioning(partitioning, GetNE());
4538 
4539  ElementToElementTable();
4540 
4541  FindPartitioningComponents(*el_to_el, _partitioning, component, num_comp);
4542 
4543  n_empty = n_mcomp = 0;
4544  for (i = 0; i < num_comp.Size(); i++)
4545  if (num_comp[i] == 0)
4546  n_empty++;
4547  else if (num_comp[i] > 1)
4548  n_mcomp++;
4549 
4550  if (n_empty > 0)
4551  {
4552  cout << "Mesh::CheckPartitioning(...) :\n"
4553  << "The following subdomains are empty :\n";
4554  for (i = 0; i < num_comp.Size(); i++)
4555  if (num_comp[i] == 0)
4556  cout << ' ' << i;
4557  cout << endl;
4558  }
4559  if (n_mcomp > 0)
4560  {
4561  cout << "Mesh::CheckPartitioning(...) :\n"
4562  << "The following subdomains are NOT connected :\n";
4563  for (i = 0; i < num_comp.Size(); i++)
4564  if (num_comp[i] > 1)
4565  cout << ' ' << i;
4566  cout << endl;
4567  }
4568  if (n_empty == 0 && n_mcomp == 0)
4569  cout << "Mesh::CheckPartitioning(...) : "
4570  "All subdomains are connected." << endl;
4571 
4572  if (el_to_el)
4573  delete el_to_el;
4574  el_to_el = NULL;
4575 }
4576 
4577 // compute the coefficients of the polynomial in t:
4578 // c(0)+c(1)*t+...+c(d)*t^d = det(A+t*B)
4579 // where A, B are (d x d), d=2,3
4580 void DetOfLinComb(const DenseMatrix &A, const DenseMatrix &B, Vector &c)
4581 {
4582  const double *a = A.Data();
4583  const double *b = B.Data();
4584 
4585  c.SetSize(A.Width()+1);
4586  switch (A.Width())
4587  {
4588  case 2:
4589  {
4590  // det(A+t*B) = |a0 a2| / |a0 b2| + |b0 a2| \ |b0 b2|
4591  // |a1 a3| + \ |a1 b3| |b1 a3| / * t + |b1 b3| * t^2
4592  c(0) = a[0]*a[3]-a[1]*a[2];
4593  c(1) = a[0]*b[3]-a[1]*b[2]+b[0]*a[3]-b[1]*a[2];
4594  c(2) = b[0]*b[3]-b[1]*b[2];
4595  }
4596  break;
4597 
4598  case 3:
4599  {
4600  /* |a0 a3 a6|
4601  * det(A+t*B) = |a1 a4 a7| +
4602  * |a2 a5 a8|
4603 
4604  * / |b0 a3 a6| |a0 b3 a6| |a0 a3 b6| \
4605  * + | |b1 a4 a7| + |a1 b4 a7| + |a1 a4 b7| | * t +
4606  * \ |b2 a5 a8| |a2 b5 a8| |a2 a5 b8| /
4607 
4608  * / |a0 b3 b6| |b0 a3 b6| |b0 b3 a6| \
4609  * + | |a1 b4 b7| + |b1 a4 b7| + |b1 b4 a7| | * t^2 +
4610  * \ |a2 b5 b8| |b2 a5 b8| |b2 b5 a8| /
4611 
4612  * |b0 b3 b6|
4613  * + |b1 b4 b7| * t^3
4614  * |b2 b5 b8| */
4615  c(0) = (a[0] * (a[4] * a[8] - a[5] * a[7]) +
4616  a[1] * (a[5] * a[6] - a[3] * a[8]) +
4617  a[2] * (a[3] * a[7] - a[4] * a[6]));
4618 
4619  c(1) = (b[0] * (a[4] * a[8] - a[5] * a[7]) +
4620  b[1] * (a[5] * a[6] - a[3] * a[8]) +
4621  b[2] * (a[3] * a[7] - a[4] * a[6]) +
4622 
4623  a[0] * (b[4] * a[8] - b[5] * a[7]) +
4624  a[1] * (b[5] * a[6] - b[3] * a[8]) +
4625  a[2] * (b[3] * a[7] - b[4] * a[6]) +
4626 
4627  a[0] * (a[4] * b[8] - a[5] * b[7]) +
4628  a[1] * (a[5] * b[6] - a[3] * b[8]) +
4629  a[2] * (a[3] * b[7] - a[4] * b[6]));
4630 
4631  c(2) = (a[0] * (b[4] * b[8] - b[5] * b[7]) +
4632  a[1] * (b[5] * b[6] - b[3] * b[8]) +
4633  a[2] * (b[3] * b[7] - b[4] * b[6]) +
4634 
4635  b[0] * (a[4] * b[8] - a[5] * b[7]) +
4636  b[1] * (a[5] * b[6] - a[3] * b[8]) +
4637  b[2] * (a[3] * b[7] - a[4] * b[6]) +
4638 
4639  b[0] * (b[4] * a[8] - b[5] * a[7]) +
4640  b[1] * (b[5] * a[6] - b[3] * a[8]) +
4641  b[2] * (b[3] * a[7] - b[4] * a[6]));
4642 
4643  c(3) = (b[0] * (b[4] * b[8] - b[5] * b[7]) +
4644  b[1] * (b[5] * b[6] - b[3] * b[8]) +
4645  b[2] * (b[3] * b[7] - b[4] * b[6]));
4646  }
4647  break;
4648 
4649  default:
4650  mfem_error("DetOfLinComb(...)");
4651  }
4652 }
4653 
4654 // compute the real roots of
4655 // z(0)+z(1)*x+...+z(d)*x^d = 0, d=2,3;
4656 // the roots are returned in x, sorted in increasing order;
4657 // it is assumed that x is at least of size d;
4658 // return the number of roots counting multiplicity;
4659 // return -1 if all z(i) are 0.
4660 int FindRoots(const Vector &z, Vector &x)
4661 {
4662  int d = z.Size()-1;
4663  if (d > 3 || d < 0)
4664  mfem_error("FindRoots(...)");
4665 
4666  while (z(d) == 0.0)
4667  {
4668  if (d == 0)
4669  return(-1);
4670  d--;
4671  }
4672  switch (d)
4673  {
4674  case 0:
4675  {
4676  return 0;
4677  }
4678 
4679  case 1:
4680  {
4681  x(0) = -z(0)/z(1);
4682  return 1;
4683  }
4684 
4685  case 2:
4686  {
4687  double a = z(2), b = z(1), c = z(0);
4688  double D = b*b-4*a*c;
4689  if (D < 0.0)
4690  {
4691  return 0;
4692  }
4693  if (D == 0.0)
4694  {
4695  x(0) = x(1) = -0.5 * b / a;
4696  return 2; // root with multiplicity 2
4697  }
4698  if (b == 0.0)
4699  {
4700  x(0) = -(x(1) = fabs(0.5 * sqrt(D) / a));
4701  return 2;
4702  }
4703  else
4704  {
4705  double t;
4706  if (b > 0.0)
4707  t = -0.5 * (b + sqrt(D));
4708  else
4709  t = -0.5 * (b - sqrt(D));
4710  x(0) = t / a;
4711  x(1) = c / t;
4712  if (x(0) > x(1))
4713  Swap<double>(x(0), x(1));
4714  return 2;
4715  }
4716  }
4717 
4718  case 3:
4719  {
4720  double a = z(2)/z(3), b = z(1)/z(3), c = z(0)/z(3);
4721 
4722  // find the real roots of x^3 + a x^2 + b x + c = 0
4723  double Q = (a * a - 3 * b) / 9;
4724  double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
4725  double Q3 = Q * Q * Q;
4726  double R2 = R * R;
4727 
4728  if (R2 == Q3)
4729  {
4730  if (Q == 0)
4731  {
4732  x(0) = x(1) = x(2) = - a / 3;
4733  }
4734  else
4735  {
4736  double sqrtQ = sqrt(Q);
4737 
4738  if (R > 0)
4739  {
4740  x(0) = -2 * sqrtQ - a / 3;
4741  x(1) = x(2) = sqrtQ - a / 3;
4742  }
4743  else
4744  {
4745  x(0) = x(1) = - sqrtQ - a / 3;
4746  x(2) = 2 * sqrtQ - a / 3;
4747  }
4748  }
4749  return 3;
4750  }
4751  else if (R2 < Q3)
4752  {
4753  double theta = acos(R / sqrt(Q3));
4754  double A = -2 * sqrt(Q);
4755  double x0, x1, x2;
4756  x0 = A * cos(theta / 3) - a / 3;
4757  x1 = A * cos((theta + 2.0 * M_PI) / 3) - a / 3;
4758  x2 = A * cos((theta - 2.0 * M_PI) / 3) - a / 3;
4759 
4760  /* Sort x0, x1, x2 */
4761  if (x0 > x1)
4762  Swap<double>(x0, x1);
4763  if (x1 > x2)
4764  {
4765  Swap<double>(x1, x2);
4766  if (x0 > x1)
4767  Swap<double>(x0, x1);
4768  }
4769  x(0) = x0;
4770  x(1) = x1;
4771  x(2) = x2;
4772  return 3;
4773  }
4774  else
4775  {
4776  double A;
4777  if (R >= 0.0)
4778  A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
4779  else
4780  A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
4781  x(0) = A + Q / A - a / 3;
4782  return 1;
4783  }
4784  }
4785  }
4786  return 0;
4787 }
4788 
4789 void FindTMax(Vector &c, Vector &x, double &tmax,
4790  const double factor, const int Dim)
4791 {
4792  const double c0 = c(0);
4793  c(0) = c0 * (1.0 - pow(factor, -Dim));
4794  int nr = FindRoots(c, x);
4795  for (int j = 0; j < nr; j++)
4796  {
4797  if (x(j) > tmax)
4798  break;
4799  if (x(j) >= 0.0)
4800  {
4801  tmax = x(j);
4802  break;
4803  }
4804  }
4805  c(0) = c0 * (1.0 - pow(factor, Dim));
4806  nr = FindRoots(c, x);
4807  for (int j = 0; j < nr; j++)
4808  {
4809  if (x(j) > tmax)
4810  break;
4811  if (x(j) >= 0.0)
4812  {
4813  tmax = x(j);
4814  break;
4815  }
4816  }
4817 }
4818 
4819 void Mesh::CheckDisplacements(const Vector &displacements, double &tmax)
4820 {
4821  int nvs = vertices.Size();
4822  DenseMatrix P, V, DS, PDS(spaceDim), VDS(spaceDim);
4823  Vector c(spaceDim+1), x(spaceDim);
4824  const double factor = 2.0;
4825 
4826  // check for tangling assuming constant speed
4827  if (tmax < 1.0)
4828  tmax = 1.0;
4829  for (int i = 0; i < NumOfElements; i++)
4830  {
4831  Element *el = elements[i];
4832  int nv = el->GetNVertices();
4833  int *v = el->GetVertices();
4834  P.SetSize(spaceDim, nv);
4835  V.SetSize(spaceDim, nv);
4836  for (int j = 0; j < spaceDim; j++)
4837  for (int k = 0; k < nv; k++)
4838  {
4839  P(j, k) = vertices[v[k]](j);
4840  V(j, k) = displacements(v[k]+j*nvs);
4841  }
4842  DS.SetSize(nv, spaceDim);
4843  const FiniteElement *fe =
4844  GetTransformationFEforElementType(el->GetType());
4845  // check if det(P.DShape+t*V.DShape) > 0 for all x and 0<=t<=1
4846  switch (el->GetType())
4847  {
4848  case Element::TRIANGLE:
4849  case Element::TETRAHEDRON:
4850  {
4851  // DS is constant
4852  fe->CalcDShape(Geometries.GetCenter(fe->GetGeomType()), DS);
4853  Mult(P, DS, PDS);
4854  Mult(V, DS, VDS);
4855  DetOfLinComb(PDS, VDS, c);
4856  if (c(0) <= 0.0)
4857  tmax = 0.0;
4858  else
4859  FindTMax(c, x, tmax, factor, Dim);
4860  }
4861  break;
4862 
4864  {
4865  const IntegrationRule &ir = fe->GetNodes();
4866  for (int j = 0; j < nv; j++)
4867  {
4868  fe->CalcDShape(ir.IntPoint(j), DS);
4869  Mult(P, DS, PDS);
4870  Mult(V, DS, VDS);
4871  DetOfLinComb(PDS, VDS, c);
4872  if (c(0) <= 0.0)
4873  tmax = 0.0;
4874  else
4875  FindTMax(c, x, tmax, factor, Dim);
4876  }
4877  }
4878  break;
4879 
4880  default:
4881  mfem_error("Mesh::CheckDisplacements(...)");
4882  }
4883  }
4884 }
4885 
4886 void Mesh::MoveVertices(const Vector &displacements)
4887 {
4888  for (int i = 0, nv = vertices.Size(); i < nv; i++)
4889  for (int j = 0; j < spaceDim; j++)
4890  vertices[i](j) += displacements(j*nv+i);
4891 }
4892 
4893 void Mesh::GetVertices(Vector &vert_coord) const
4894 {
4895  int nv = vertices.Size();
4896  vert_coord.SetSize(nv*spaceDim);
4897  for (int i = 0; i < nv; i++)
4898  for (int j = 0; j < spaceDim; j++)
4899  vert_coord(j*nv+i) = vertices[i](j);
4900 }
4901 
4902 void Mesh::SetVertices(const Vector &vert_coord)
4903 {
4904  for (int i = 0, nv = vertices.Size(); i < nv; i++)
4905  for (int j = 0; j < spaceDim; j++)
4906  vertices[i](j) = vert_coord(j*nv+i);
4907 }
4908 
4909 void Mesh::GetNode(int i, double *coord)
4910 {
4911  if (Nodes)
4912  {
4913  FiniteElementSpace *fes = Nodes->FESpace();
4914  for (int j = 0; j < spaceDim; j++)
4915  coord[j] = (*Nodes)(fes->DofToVDof(i, j));
4916  }
4917  else
4918  {
4919  for (int j = 0; j < spaceDim; j++)
4920  coord[j] = vertices[i](j);
4921  }
4922 }
4923 
4924 void Mesh::SetNode(int i, const double *coord)
4925 {
4926  if (Nodes)
4927  {
4928  FiniteElementSpace *fes = Nodes->FESpace();
4929  for (int j = 0; j < spaceDim; j++)
4930  (*Nodes)(fes->DofToVDof(i, j)) = coord[j];
4931  }
4932  else
4933  {
4934  for (int j = 0; j < spaceDim; j++)
4935  vertices[i](j) = coord[j];
4936 
4937  }
4938 }
4939 
4940 void Mesh::MoveNodes(const Vector &displacements)
4941 {
4942  if (Nodes)
4943  (*Nodes) += displacements;
4944  else
4945  MoveVertices(displacements);
4946 }
4947 
4948 void Mesh::GetNodes(Vector &node_coord) const
4949 {
4950  if (Nodes)
4951  node_coord = (*Nodes);
4952  else
4953  GetVertices(node_coord);
4954 }
4955 
4956 void Mesh::SetNodes(const Vector &node_coord)
4957 {
4958  if (Nodes)
4959  (*Nodes) = node_coord;
4960  else
4961  SetVertices(node_coord);
4962 }
4963 
4964 void Mesh::NewNodes(GridFunction &nodes, bool make_owner)
4965 {
4966  if (own_nodes) delete Nodes;
4967  Nodes = &nodes;
4968  spaceDim = Nodes->FESpace()->GetVDim();
4969  own_nodes = (int)make_owner;
4970 
4971  if (NURBSext != nodes.FESpace()->GetNURBSext())
4972  {
4973  delete NURBSext;
4974  NURBSext = nodes.FESpace()->StealNURBSext();
4975  }
4976 }
4977 
4978 void Mesh::SwapNodes(GridFunction *&nodes, int &own_nodes_)
4979 {
4980  mfem::Swap<GridFunction*>(Nodes, nodes);
4981  mfem::Swap<int>(own_nodes, own_nodes_);
4982  // TODO:
4983  // if (nodes)
4984  // nodes->FESpace()->MakeNURBSextOwner();
4985  // NURBSext = (Nodes) ? Nodes->FESpace()->StealNURBSext() : NULL;
4986 }
4987 
4988 void Mesh::AverageVertices(int * indexes, int n, int result)
4989 {
4990  int j, k;
4991 
4992  for (k = 0; k < spaceDim; k++)
4993  vertices[result](k) = vertices[indexes[0]](k);
4994 
4995  for (j = 1; j < n; j++)
4996  for (k = 0; k < spaceDim; k++)
4997  vertices[result](k) += vertices[indexes[j]](k);
4998 
4999  for (k = 0; k < spaceDim; k++)
5000  vertices[result](k) *= (1.0 / n);
5001 }
5002 
5004 {
5005  Nodes->FESpace()->UpdateAndInterpolate(Nodes);
5006  // update the vertices?
5007 }
5008 
5010 {
5011  int i, j, *v, vv[2], attr, wtls = WantTwoLevelState;
5012  const int *e;
5013 
5014  if (Nodes) // curved mesh
5015  {
5016  UseTwoLevelState(1);
5017  }
5018 
5019  SetState(Mesh::NORMAL);
5020 
5021  if (el_to_edge == NULL)
5022  {
5023  el_to_edge = new Table;
5024  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5025  }
5026 
5027  int oedge = NumOfVertices;
5028  int oelem = oedge + NumOfEdges;
5029 
5030  DeleteCoarseTables();
5031 
5032  vertices.SetSize(oelem + NumOfElements);
5033 
5034  for (i = 0; i < NumOfElements; i++)
5035  {
5036  v = elements[i]->GetVertices();
5037 
5038  AverageVertices(v, 4, oelem+i);
5039 
5040  e = el_to_edge->GetRow(i);
5041 
5042  vv[0] = v[0], vv[1] = v[1]; AverageVertices(vv, 2, oedge+e[0]);
5043  vv[0] = v[1], vv[1] = v[2]; AverageVertices(vv, 2, oedge+e[1]);
5044  vv[0] = v[2], vv[1] = v[3]; AverageVertices(vv, 2, oedge+e[2]);
5045  vv[0] = v[3], vv[1] = v[0]; AverageVertices(vv, 2, oedge+e[3]);
5046  }
5047 
5048  elements.SetSize(4 * NumOfElements);
5049  for (i = 0; i < NumOfElements; i++)
5050  {
5051  attr = elements[i]->GetAttribute();
5052  v = elements[i]->GetVertices();
5053  e = el_to_edge->GetRow(i);
5054  j = NumOfElements + 3 * i;
5055 
5056  elements[j+0] = new Quadrilateral(oedge+e[0], v[1], oedge+e[1],
5057  oelem+i, attr);
5058  elements[j+1] = new Quadrilateral(oelem+i, oedge+e[1], v[2],
5059  oedge+e[2], attr);
5060  elements[j+2] = new Quadrilateral(oedge+e[3], oelem+i, oedge+e[2],
5061  v[3], attr);
5062 
5063  if (WantTwoLevelState)
5064  {
5065  QuadrisectedElement *qe;
5066 
5067  qe = new QuadrisectedElement(elements[i]->Duplicate(this));
5068  qe->FirstChild = elements[i];
5069  qe->Child2 = j;
5070  qe->Child3 = j+1;
5071  qe->Child4 = j+2;
5072  elements[i] = qe;
5073  }
5074 
5075  v[1] = oedge+e[0];
5076  v[2] = oelem+i;
5077  v[3] = oedge+e[3];
5078  }
5079 
5080  boundary.SetSize(2 * NumOfBdrElements);
5081  for (i = 0; i < NumOfBdrElements; i++)
5082  {
5083  attr = boundary[i]->GetAttribute();
5084  v = boundary[i]->GetVertices();
5085  j = NumOfBdrElements + i;
5086 
5087  boundary[j] = new Segment(oedge+be_to_edge[i], v[1], attr);
5088 
5089  if (WantTwoLevelState)
5090  {
5091 #ifdef MFEM_USE_MEMALLOC
5092  BisectedElement *be = BEMemory.Alloc();
5093 #else
5094  BisectedElement *be = new BisectedElement;
5095 #endif
5096  be->SetCoarseElem(boundary[i]->Duplicate(this));
5097  be->FirstChild = boundary[i];
5098  be->SecondChild = j;
5099  boundary[i] = be;
5100  }
5101 
5102  v[1] = oedge+be_to_edge[i];
5103  }
5104 
5105  if (WantTwoLevelState)
5106  {
5107  c_NumOfVertices = NumOfVertices;
5108  c_NumOfEdges = NumOfEdges;
5109  c_NumOfElements = NumOfElements;
5110  c_NumOfBdrElements = NumOfBdrElements;
5111 
5113  State = Mesh::TWO_LEVEL_FINE;
5114  }
5115 
5116  NumOfVertices = oelem + NumOfElements;
5117  NumOfElements = 4 * NumOfElements;
5118  NumOfBdrElements = 2 * NumOfBdrElements;
5119  NumOfFaces = 0;
5120 
5121  if (WantTwoLevelState)
5122  {
5123  f_NumOfVertices = NumOfVertices;
5124  f_NumOfElements = NumOfElements;
5125  f_NumOfBdrElements = NumOfBdrElements;
5126  }
5127 
5128  if (el_to_edge != NULL)
5129  {
5130  if (WantTwoLevelState)
5131  {
5132  c_el_to_edge = el_to_edge;
5133  mfem::Swap(be_to_edge, fc_be_to_edge); // save coarse be_to_edge
5134  f_el_to_edge = new Table;
5135  NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
5136  el_to_edge = f_el_to_edge;
5137  f_NumOfEdges = NumOfEdges;
5138  }
5139  else
5140  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5141  GenerateFaces();
5142  }
5143 
5144  if (Nodes) // curved mesh
5145  {
5146  UpdateNodes();
5147  UseTwoLevelState(wtls);
5148  }
5149 
5150 #ifdef MFEM_DEBUG
5151  CheckElementOrientation(false);
5152  CheckBdrElementOrientation(false);
5153 #endif
5154 }
5155 
5157 {
5158  int i, wtls = WantTwoLevelState;
5159  int * v;
5160  const int *e, *f;
5161  int vv[4];
5162 
5163  if (Nodes) // curved mesh
5164  {
5165  UseTwoLevelState(1);
5166  }
5167 
5168  SetState(Mesh::NORMAL);
5169 
5170  if (el_to_edge == NULL)
5171  {
5172  el_to_edge = new Table;
5173  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5174  }
5175  if (el_to_face == NULL)
5176  GetElementToFaceTable();
5177 
5178  int oedge = NumOfVertices;
5179  int oface = oedge + NumOfEdges;
5180  int oelem = oface + NumOfFaces;
5181 
5182  DeleteCoarseTables();
5183 
5184  vertices.SetSize(oelem + NumOfElements);
5185  for (i = 0; i < NumOfElements; i++)
5186  {
5187  MFEM_ASSERT(elements[i]->GetType() == Element::HEXAHEDRON,
5188  "Element is not a hex!");
5189  v = elements[i]->GetVertices();
5190 
5191  AverageVertices(v, 8, oelem+i);
5192 
5193  f = el_to_face->GetRow(i);
5194 
5195  for (int j = 0; j < 6; j++)
5196  {
5197  for (int k = 0; k < 4; k++)
5198  vv[k] = v[hex_faces[j][k]];
5199  AverageVertices(vv, 4, oface+f[j]);
5200  }
5201 
5202  e = el_to_edge->GetRow(i);
5203 
5204  for (int j = 0; j < 12; j++)
5205  {
5206  for (int k = 0; k < 2; k++)
5207  vv[k] = v[Hexahedron::edges[j][k]];
5208  AverageVertices(vv, 2, oedge+e[j]);
5209  }
5210  }
5211 
5212  int attr, j, k;
5213  elements.SetSize(8 * NumOfElements);
5214  for (i = 0; i < NumOfElements; i++)
5215  {
5216  attr = elements[i]->GetAttribute();
5217  v = elements[i]->GetVertices();
5218  e = el_to_edge->GetRow(i);
5219  f = el_to_face->GetRow(i);
5220  j = NumOfElements + 7 * i;
5221 
5222  elements[j+0] = new Hexahedron(oedge+e[0], v[1], oedge+e[1], oface+f[0],
5223  oface+f[1], oedge+e[9], oface+f[2],
5224  oelem+i, attr);
5225  elements[j+1] = new Hexahedron(oface+f[0], oedge+e[1], v[2], oedge+e[2],
5226  oelem+i, oface+f[2], oedge+e[10],
5227  oface+f[3], attr);
5228  elements[j+2] = new Hexahedron(oedge+e[3], oface+f[0], oedge+e[2], v[3],
5229  oface+f[4], oelem+i, oface+f[3],
5230  oedge+e[11], attr);
5231  elements[j+3] = new Hexahedron(oedge+e[8], oface+f[1], oelem+i,
5232  oface+f[4], v[4], oedge+e[4], oface+f[5],
5233  oedge+e[7], attr);
5234  elements[j+4] = new Hexahedron(oface+f[1], oedge+e[9], oface+f[2],
5235  oelem+i, oedge+e[4], v[5], oedge+e[5],
5236  oface+f[5], attr);
5237  elements[j+5] = new Hexahedron(oelem+i, oface+f[2], oedge+e[10],
5238  oface+f[3], oface+f[5], oedge+e[5], v[6],
5239  oedge+e[6], attr);
5240  elements[j+6] = new Hexahedron(oface+f[4], oelem+i, oface+f[3],
5241  oedge+e[11], oedge+e[7], oface+f[5],
5242  oedge+e[6], v[7], attr);
5243 
5244  if (WantTwoLevelState)
5245  {
5246  OctasectedElement *oe;
5247 
5248  oe = new OctasectedElement(elements[i]->Duplicate(this));
5249  oe->FirstChild = elements[i];
5250  for (k = 0; k < 7; k++)
5251  oe->Child[k] = j + k;
5252  elements[i] = oe;
5253  }
5254 
5255  v[1] = oedge+e[0];
5256  v[2] = oface+f[0];
5257  v[3] = oedge+e[3];
5258  v[4] = oedge+e[8];
5259  v[5] = oface+f[1];
5260  v[6] = oelem+i;
5261  v[7] = oface+f[4];
5262  }
5263 
5264  boundary.SetSize(4 * NumOfBdrElements);
5265  for (i = 0; i < NumOfBdrElements; i++)
5266  {
5267  MFEM_ASSERT(boundary[i]->GetType() == Element::QUADRILATERAL,
5268  "boundary Element is not a quad!");
5269  attr = boundary[i]->GetAttribute();
5270  v = boundary[i]->GetVertices();
5271  e = bel_to_edge->GetRow(i);
5272  f = & be_to_face[i];
5273  j = NumOfBdrElements + 3 * i;
5274 
5275  boundary[j+0] = new Quadrilateral(oedge+e[0], v[1], oedge+e[1],
5276  oface+f[0], attr);
5277  boundary[j+1] = new Quadrilateral(oface+f[0], oedge+e[1], v[2],
5278  oedge+e[2], attr);
5279  boundary[j+2] = new Quadrilateral(oedge+e[3], oface+f[0], oedge+e[2],
5280  v[3], attr);
5281 
5282  if (WantTwoLevelState)
5283  {
5284  QuadrisectedElement *qe;
5285 
5286  qe = new QuadrisectedElement(boundary[i]->Duplicate(this));
5287  qe->FirstChild = boundary[i];
5288  qe->Child2 = j;
5289  qe->Child3 = j+1;
5290  qe->Child4 = j+2;
5291  boundary[i] = qe;
5292  }
5293 
5294  v[1] = oedge+e[0];
5295  v[2] = oface+f[0];
5296  v[3] = oedge+e[3];
5297  }
5298 
5299  if (WantTwoLevelState)
5300  {
5301  c_NumOfVertices = NumOfVertices;
5302  c_NumOfEdges = NumOfEdges;
5303  c_NumOfFaces = NumOfFaces;
5304  c_NumOfElements = NumOfElements;
5305  c_NumOfBdrElements = NumOfBdrElements;
5306 
5308  State = Mesh::TWO_LEVEL_FINE;
5309  }
5310 
5311  NumOfVertices = oelem + NumOfElements;
5312  NumOfElements = 8 * NumOfElements;
5313  NumOfBdrElements = 4 * NumOfBdrElements;
5314 
5315  if (WantTwoLevelState)
5316  {
5317  f_NumOfVertices = NumOfVertices;
5318  f_NumOfElements = NumOfElements;
5319  f_NumOfBdrElements = NumOfBdrElements;
5320  }
5321 
5322  if (el_to_face != NULL)
5323  {
5324  if (WantTwoLevelState)
5325  {
5326  c_el_to_face = el_to_face;
5327  el_to_face = NULL;
5328  mfem::Swap(faces_info, fc_faces_info);
5329  }
5330  GetElementToFaceTable();
5331  GenerateFaces();
5332  if (WantTwoLevelState)
5333  {
5334  f_el_to_face = el_to_face;
5335  f_NumOfFaces = NumOfFaces;
5336  }
5337  }
5338 
5339 #ifdef MFEM_DEBUG
5340  CheckBdrElementOrientation(false);
5341 #endif
5342 
5343  if (el_to_edge != NULL)
5344  {
5345  if (WantTwoLevelState)
5346  {
5347  c_el_to_edge = el_to_edge;
5348  f_el_to_edge = new Table;
5349  c_bel_to_edge = bel_to_edge;
5350  bel_to_edge = NULL;
5351  NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
5352  el_to_edge = f_el_to_edge;
5353  f_bel_to_edge = bel_to_edge;
5354  f_NumOfEdges = NumOfEdges;
5355  }
5356  else
5357  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5358  }
5359 
5360  if (Nodes) // curved mesh
5361  {
5362  UpdateNodes();
5363  UseTwoLevelState(wtls);
5364  }
5365 
5366  // When 'WantTwoLevelState' is true the coarse level
5367  // 'be_to_face' and 'faces'
5368  // are destroyed !!!
5369 }
5370 
5371 void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
5372 {
5373  int i, j, ind, nedges, wtls = WantTwoLevelState;
5374  Array<int> v;
5375 
5376  if (ncmesh)
5377  {
5378  mfem_error("Local and nonconforming refinements cannot be mixed.");
5379  }
5380 
5381  if (Nodes) // curved mesh
5382  {
5383  UseTwoLevelState(1);
5384  }
5385 
5386  SetState(Mesh::NORMAL);
5387  DeleteCoarseTables();
5388 
5389  if (Dim == 1) // --------------------------------------------------------
5390  {
5391  if (WantTwoLevelState)
5392  {
5393  c_NumOfVertices = NumOfVertices;
5394  c_NumOfElements = NumOfElements;
5395  c_NumOfBdrElements = NumOfBdrElements;
5396  c_NumOfEdges = 0;
5397  }
5398  int cne = NumOfElements, cnv = NumOfVertices;
5399  NumOfVertices += marked_el.Size();
5400  NumOfElements += marked_el.Size();
5401  vertices.SetSize(NumOfVertices);
5402  elements.SetSize(NumOfElements);
5403  for (j = 0; j < marked_el.Size(); j++)
5404  {
5405  i = marked_el[j];
5406  Segment *c_seg = (Segment *)elements[i];
5407  int *vert = c_seg->GetVertices(), attr = c_seg->GetAttribute();
5408  int new_v = cnv + j, new_e = cne + j;
5409  AverageVertices(vert, 2, new_v);
5410  elements[new_e] = new Segment(new_v, vert[1], attr);
5411  if (WantTwoLevelState)
5412  {
5413 #ifdef MFEM_USE_MEMALLOC
5414  BisectedElement *aux = BEMemory.Alloc();
5415  aux->SetCoarseElem(c_seg);
5416 #else
5417  BisectedElement *aux = new BisectedElement(c_seg);
5418 #endif
5419  aux->FirstChild = new Segment(vert[0], new_v, attr);
5420  aux->SecondChild = new_e;
5421  elements[i] = aux;
5422  }
5423  else
5424  {
5425  vert[1] = new_v;
5426  }
5427  }
5428  if (WantTwoLevelState)
5429  {
5430  f_NumOfVertices = NumOfVertices;
5431  f_NumOfElements = NumOfElements;
5432  f_NumOfBdrElements = NumOfBdrElements;
5433  f_NumOfEdges = 0;
5434 
5436  State = Mesh::TWO_LEVEL_FINE;
5437  }
5438  GenerateFaces();
5439  } // end of 'if (Dim == 1)'
5440  else if (Dim == 2) // ---------------------------------------------------
5441  {
5442  if (WantTwoLevelState)
5443  {
5444  c_NumOfVertices = NumOfVertices;
5445  c_NumOfEdges = NumOfEdges;
5446  c_NumOfElements = NumOfElements;
5447  c_NumOfBdrElements = NumOfBdrElements;
5448  }
5449 
5450  // 1. Get table of vertex to vertex connections.
5451  DSTable v_to_v(NumOfVertices);
5452  GetVertexToVertexTable(v_to_v);
5453 
5454  // 2. Get edge to element connections in arrays edge1 and edge2
5455  nedges = v_to_v.NumberOfEntries();
5456  int *edge1 = new int[nedges];
5457  int *edge2 = new int[nedges];
5458  int *middle = new int[nedges];
5459 
5460  for (i = 0; i < nedges; i++)
5461  edge1[i] = edge2[i] = middle[i] = -1;
5462 
5463  for (i = 0; i < NumOfElements; i++)
5464  {
5465  elements[i]->GetVertices(v);
5466  for (j = 1; j < v.Size(); j++)
5467  {
5468  ind = v_to_v(v[j-1], v[j]);
5469  (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
5470  }
5471  ind = v_to_v(v[0], v[v.Size()-1]);
5472  (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
5473  }
5474 
5475  // 3. Do the red refinement.
5476  for (i = 0; i < marked_el.Size(); i++)
5477  RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
5478 
5479  // 4. Do the green refinement (to get conforming mesh).
5480  int need_refinement;
5481  do
5482  {
5483  need_refinement = 0;
5484  for (i = 0; i < nedges; i++)
5485  if (middle[i] != -1 && edge1[i] != -1)
5486  {
5487  need_refinement = 1;
5488  GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
5489  }
5490  }
5491  while (need_refinement == 1);
5492 
5493  // 5. Update the boundary elements.
5494  int v1[2], v2[2], bisect, temp;
5495  temp = NumOfBdrElements;
5496  for (i = 0; i < temp; i++)
5497  {
5498  boundary[i]->GetVertices(v);
5499  bisect = v_to_v(v[0], v[1]);
5500  if (middle[bisect] != -1) // the element was refined (needs updating)
5501  {
5502  if (boundary[i]->GetType() == Element::SEGMENT)
5503  {
5504  v1[0] = v[0]; v1[1] = middle[bisect];
5505  v2[0] = middle[bisect]; v2[1] = v[1];
5506 
5507  if (WantTwoLevelState)
5508  {
5509  boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
5510 #ifdef MFEM_USE_MEMALLOC
5511  BisectedElement *aux = BEMemory.Alloc();
5512  aux->SetCoarseElem(boundary[i]);
5513 #else
5514  BisectedElement *aux = new BisectedElement(boundary[i]);
5515 #endif
5516  aux->FirstChild =
5517  new Segment(v1, boundary[i]->GetAttribute());
5518  aux->SecondChild = NumOfBdrElements;
5519  boundary[i] = aux;
5520  NumOfBdrElements++;
5521  }
5522  else
5523  {
5524  boundary[i]->SetVertices(v1);
5525  boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
5526  }
5527  }
5528  else
5529  mfem_error("Only bisection of segment is implemented"
5530  " for bdr elem.");
5531  }
5532  }
5533  NumOfBdrElements = boundary.Size();
5534 
5535  // 6. Free the allocated memory.
5536  delete [] edge1;
5537  delete [] edge2;
5538  delete [] middle;
5539 
5540  if (WantTwoLevelState)
5541  {
5542  f_NumOfVertices = NumOfVertices;
5543  f_NumOfElements = NumOfElements;
5544  f_NumOfBdrElements = NumOfBdrElements;
5546  State = Mesh::TWO_LEVEL_FINE;
5547  }
5548 
5549  if (el_to_edge != NULL)
5550  {
5551  if (WantTwoLevelState)
5552  {
5553  c_el_to_edge = el_to_edge;
5554  mfem::Swap(be_to_edge, fc_be_to_edge); // save coarse be_to_edge
5555  f_el_to_edge = new Table;
5556  NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
5557  el_to_edge = f_el_to_edge;
5558  f_NumOfEdges = NumOfEdges;
5559  }
5560  else
5561  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5562  GenerateFaces();
5563  }
5564 
5565  }
5566  else if (Dim == 3) // ---------------------------------------------------
5567  {
5568  if (WantTwoLevelState)
5569  {
5570  c_NumOfVertices = NumOfVertices;
5571  c_NumOfEdges = NumOfEdges;
5572  c_NumOfFaces = NumOfFaces;
5573  c_NumOfElements = NumOfElements;
5574  c_NumOfBdrElements = NumOfBdrElements;
5575  }
5576 
5577  // 1. Get table of vertex to vertex connections.
5578  DSTable v_to_v(NumOfVertices);
5579  GetVertexToVertexTable(v_to_v);
5580 
5581  // 2. Get edge to element connections in arrays edge1 and edge2
5582  nedges = v_to_v.NumberOfEntries();
5583  int *middle = new int[nedges];
5584 
5585  for (i = 0; i < nedges; i++)
5586  middle[i] = -1;
5587 
5588  // 3. Do the red refinement.
5589  int ii;
5590  switch (type)
5591  {
5592  case 1:
5593  for (i = 0; i < marked_el.Size(); i++)
5594  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5595  break;
5596  case 2:
5597  for (i = 0; i < marked_el.Size(); i++)
5598  {
5599  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5600 
5601  Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
5602  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5603  }
5604  break;
5605  case 3:
5606  for (i = 0; i < marked_el.Size(); i++)
5607  {
5608  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5609 
5610  ii = NumOfElements - 1;
5611  Bisection(ii, v_to_v, NULL, NULL, middle);
5612  Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
5613  Bisection(ii, v_to_v, NULL, NULL, middle);
5614 
5615  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5616  Bisection(NumOfElements-1, v_to_v, NULL, NULL, middle);
5617  Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5618  }
5619  break;
5620  }
5621 
5622  if (WantTwoLevelState)
5623  {
5625  State = Mesh::TWO_LEVEL_FINE;
5626  }
5627 
5628  // 4. Do the green refinement (to get conforming mesh).
5629  int need_refinement;
5630  // int need_refinement, onoe, max_gen = 0;
5631  do
5632  {
5633  // int redges[2], type, flag;
5634  need_refinement = 0;
5635  // onoe = NumOfElements;
5636  // for (i = 0; i < onoe; i++)
5637  for (i = 0; i < NumOfElements; i++)
5638  {
5639  // ((Tetrahedron *)elements[i])->ParseRefinementFlag(redges, type, flag);
5640  // if (flag > max_gen) max_gen = flag;
5641  if (elements[i]->NeedRefinement(v_to_v, middle))
5642  {
5643  need_refinement = 1;
5644  Bisection(i, v_to_v, NULL, NULL, middle);
5645  }
5646  }
5647  }
5648  while (need_refinement == 1);
5649 
5650  // cout << "Maximum generation: " << max_gen << endl;
5651 
5652  // 5. Update the boundary elements.
5653  do
5654  {
5655  need_refinement = 0;
5656  for (i = 0; i < NumOfBdrElements; i++)
5657  if (boundary[i]->NeedRefinement(v_to_v, middle))
5658  {
5659  need_refinement = 1;
5660  Bisection(i, v_to_v, middle);
5661  }
5662  }
5663  while (need_refinement == 1);
5664 
5665  // 6. Un-mark the Pf elements.
5666  int refinement_edges[2], type, flag;
5667  for (i = 0; i < NumOfElements; i++)
5668  {
5669  Element *El = elements[i];
5670  while (El->GetType() == Element::BISECTED)
5671  El = ((BisectedElement *) El)->FirstChild;
5672  ((Tetrahedron *)El)->ParseRefinementFlag(refinement_edges, type,
5673  flag);
5674  if (type == Tetrahedron::TYPE_PF)
5675  ((Tetrahedron *)El)->CreateRefinementFlag(refinement_edges,
5677  flag);
5678  }
5679 
5680  NumOfBdrElements = boundary.Size();
5681 
5682  // 7. Free the allocated memory.
5683  delete [] middle;
5684 
5685  if (el_to_edge != NULL)
5686  {
5687  if (WantTwoLevelState)
5688  {
5689  c_el_to_edge = el_to_edge;
5690  f_el_to_edge = new Table;
5691  c_bel_to_edge = bel_to_edge;
5692  bel_to_edge = NULL;
5693  NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
5694  el_to_edge = f_el_to_edge;
5695  f_bel_to_edge = bel_to_edge;
5696  }
5697  else
5698  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5699  }
5700  if (el_to_face != NULL)
5701  {
5702  if (WantTwoLevelState)
5703  {
5704  c_el_to_face = el_to_face;
5705  el_to_face = NULL;
5706  mfem::Swap(faces_info, fc_faces_info);
5707  }
5708  GetElementToFaceTable();
5709  GenerateFaces();
5710  if (WantTwoLevelState)
5711  {
5712  f_el_to_face = el_to_face;
5713  }
5714  }
5715 
5716  if (WantTwoLevelState)
5717  {
5718  f_NumOfVertices = NumOfVertices;
5719  f_NumOfEdges = NumOfEdges;
5720  f_NumOfFaces = NumOfFaces;
5721  f_NumOfElements = NumOfElements;
5722  f_NumOfBdrElements = NumOfBdrElements;
5723  }
5724 
5725  } // end 'if (Dim == 3)'
5726 
5727  if (Nodes) // curved mesh
5728  {
5729  UpdateNodes();
5730  UseTwoLevelState(wtls);
5731  }
5732 
5733 #ifdef MFEM_DEBUG
5734  CheckElementOrientation(false);
5735 #endif
5736 }
5737 
5739  int nc_limit)
5740 {
5741  if (NURBSext)
5742  {
5743  mfem_error("Mesh::NonconformingRefinement: NURBS meshes are not supported."
5744  " Project the NURBS to Nodes first.");
5745  }
5746 
5747  int wtls = WantTwoLevelState;
5748 
5749  if (Nodes) // curved mesh
5750  {
5751  UseTwoLevelState(1);
5752  }
5753 
5754  SetState(Mesh::NORMAL);
5755  DeleteCoarseTables();
5756 
5757  if (!ncmesh)
5758  {
5759  // start tracking refinement hierarchy
5760  ncmesh = new NCMesh(this);
5761  }
5762 
5763  if (WantTwoLevelState)
5764  ncmesh->MarkCoarseLevel();
5765 
5766  // do the refinements
5767  ncmesh->Refine(refinements);
5768 
5769  if (nc_limit > 0)
5770  ncmesh->LimitNCLevel(nc_limit);
5771 
5772  // create a second mesh containing the finest elements from 'ncmesh'
5773  Mesh* mesh2 = new Mesh(*ncmesh);
5774 
5775  ncmesh->SetEdgeIndicesFromMesh(mesh2);
5776  ncmesh->SetFaceIndicesFromMesh(mesh2);
5777 
5778  // now swap the meshes, the second mesh will become the old coarse mesh
5779  // and this mesh will be the new fine mesh
5780  Swap(*mesh2, false);
5781 
5782  // retain the coarse mesh if two-level state was requested, delete otherwise
5783  if (WantTwoLevelState)
5784  {
5785  nc_coarse_level = mesh2;
5786  State = TWO_LEVEL_FINE;
5787  }
5788  else
5789  delete mesh2;
5790 
5791  if (Nodes) // curved mesh
5792  {
5793  UpdateNodes();
5794  UseTwoLevelState(wtls);
5795  }
5796 }
5797 
5799 {
5800  Dim = ncmesh.Dimension();
5801 
5802  Init();
5803  InitTables();
5804 
5805  ncmesh.GetVerticesElementsBoundary(vertices, elements, boundary);
5806 
5807  NumOfVertices = vertices.Size();
5808  NumOfElements = elements.Size();
5809  NumOfBdrElements = boundary.Size();
5810 
5811  // set the mesh type ('meshgen')
5812  SetMeshGen();
5813 
5814  NumOfEdges = NumOfFaces = 0;
5815 
5816  if (Dim > 2)
5817  {
5818  GetElementToFaceTable();
5819  GenerateFaces();
5820 #ifdef MFEM_DEBUG
5821  CheckBdrElementOrientation(false);
5822 #endif
5823  }
5824 
5825  el_to_edge = new Table;
5826  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5827  c_el_to_edge = NULL;
5828 
5829  SetAttributes();
5830 }
5831 
5832 void Mesh::Swap(Mesh& other, bool non_geometry)
5833 {
5834  mfem::Swap(Dim, other.Dim);
5835 
5836  mfem::Swap(NumOfVertices, other.NumOfVertices);
5837  mfem::Swap(NumOfElements, other.NumOfElements);
5838  mfem::Swap(NumOfBdrElements, other.NumOfBdrElements);
5839  mfem::Swap(NumOfEdges, other.NumOfEdges);
5840  mfem::Swap(NumOfFaces, other.NumOfFaces);
5841 
5842  mfem::Swap(meshgen, other.meshgen);
5843 
5844  mfem::Swap(elements, other.elements);
5845  mfem::Swap(vertices, other.vertices);
5846  mfem::Swap(boundary, other.boundary);
5847  mfem::Swap(faces, other.faces);
5848  mfem::Swap(faces_info, other.faces_info);
5849 
5850  mfem::Swap(el_to_edge, other.el_to_edge);
5851  mfem::Swap(el_to_face, other.el_to_face);
5852  mfem::Swap(el_to_el, other.el_to_el);
5853  mfem::Swap(be_to_edge, other.be_to_edge);
5854  mfem::Swap(bel_to_edge, other.bel_to_edge);
5855  mfem::Swap(be_to_face, other.be_to_face);
5856  mfem::Swap(face_edge, other.face_edge);
5857  mfem::Swap(edge_vertex, other.edge_vertex);
5858 
5859  mfem::Swap(attributes, other.attributes);
5860  mfem::Swap(bdr_attributes, other.bdr_attributes);
5861 
5862  if (non_geometry)
5863  {
5864  mfem::Swap(NURBSext, other.NURBSext);
5865  mfem::Swap(ncmesh, other.ncmesh);
5866 
5867  mfem::Swap(Nodes, other.Nodes);
5868  mfem::Swap(own_nodes, other.own_nodes);
5869  }
5870 
5871  // NOTE: two-level-state related members are ignored here
5872 }
5873 
5875 {
5876  if (NURBSext)
5877  NURBSUniformRefinement();
5878  else if (meshgen == 1)
5879  {
5880  Array<int> elem_to_refine(GetNE());
5881 
5882  for (int i = 0; i < elem_to_refine.Size(); i++)
5883  elem_to_refine[i] = i;
5884  LocalRefinement(elem_to_refine);
5885  }
5886  else if (Dim == 2)
5887  QuadUniformRefinement();
5888  else if (Dim == 3)
5889  HexUniformRefinement();
5890  else
5891  mfem_error("Mesh::UniformRefinement()");
5892 }
5893 
5894 void Mesh::GeneralRefinement(Array<Refinement> &refinements, int nonconforming,
5895  int nc_limit)
5896 {
5897  if (nonconforming < 0)
5898  {
5899  // determine if nonconforming refinement is suitable
5900  int type = elements[0]->GetType();
5901  if (type == Element::HEXAHEDRON || type == Element::QUADRILATERAL)
5902  nonconforming = 1;
5903  else
5904  nonconforming = 0;
5905  }
5906 
5907  if (nonconforming)
5908  {
5909  // non-conforming refinement (hanging nodes)
5910  NonconformingRefinement(refinements, nc_limit);
5911  }
5912  else
5913  {
5914  Array<int> el_to_refine;
5915  for (int i = 0; i < refinements.Size(); i++)
5916  el_to_refine.Append(refinements[i].index);
5917 
5918  // infer 'type' of local refinement from first element's 'ref_type'
5919  int type, rt = (refinements.Size() ? refinements[0].ref_type : 7);
5920  if (rt == 1 || rt == 2 || rt == 4)
5921  type = 1;
5922  else if (rt == 3 || rt == 5 || rt == 6)
5923  type = 2;
5924  else
5925  type = 3;
5926 
5927  // red-green refinement, no hanging nodes
5928  LocalRefinement(el_to_refine, type);
5929  }
5930 }
5931 
5932 void Mesh::GeneralRefinement(Array<int> &el_to_refine, int nonconforming,
5933  int nc_limit)
5934 {
5935  Array<Refinement> refinements;
5936  for (int i = 0; i < el_to_refine.Size(); i++)
5937  refinements.Append(Refinement(el_to_refine[i], 7));
5938 
5939  GeneralRefinement(refinements, nonconforming, nc_limit);
5940 }
5941 
5942 void Mesh::Bisection(int i, const DSTable &v_to_v,
5943  int *edge1, int *edge2, int *middle)
5944 {
5945  int *vert;
5946  int v[2][4], v_new, bisect, t;
5947  Element **pce = &elements[i];
5948  Vertex V;
5949 
5950  t = pce[0]->GetType();
5951  if (WantTwoLevelState)
5952  {
5953  while (1)
5954  {
5955  if (t == Element::BISECTED)
5956  pce = & ( ((BisectedElement *) pce[0])->FirstChild );
5957  else if (t == Element::QUADRISECTED)
5958  pce = & ( ((QuadrisectedElement *) pce[0])->FirstChild );
5959  else
5960  break;
5961  t = pce[0]->GetType();
5962  }
5963  }
5964 
5965  if (t == Element::TRIANGLE)
5966  {
5967  Triangle *tri = (Triangle *) pce[0];
5968 
5969  vert = tri->GetVertices();
5970 
5971  // 1. Get the index for the new vertex in v_new.
5972  bisect = v_to_v(vert[0], vert[1]);
5973 #ifdef MFEM_DEBUG
5974  if (bisect < 0)
5975  mfem_error("Mesh::Bisection(...) of triangle! #1");
5976 #endif
5977  if (middle[bisect] == -1)
5978  {
5979  v_new = NumOfVertices++;
5980  for (int d = 0; d < spaceDim; d++)
5981  V(d) = 0.5 * (vertices[vert[0]](d) + vertices[vert[1]](d));
5982  vertices.Append(V);
5983 
5984  // Put the element that may need refinement (because of this
5985  // bisection) in edge1, or -1 if no more refinement is needed.
5986  if (edge1[bisect] == i)
5987  edge1[bisect] = edge2[bisect];
5988 
5989  middle[bisect] = v_new;
5990  }
5991  else
5992  {
5993  v_new = middle[bisect];
5994 
5995  // This edge will require no more refinement.
5996  edge1[bisect] = -1;
5997  }
5998 
5999  // 2. Set the node indices for the new elements in v[0] and v[1] so that
6000  // the edge marked for refinement is between the first two nodes.
6001  v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
6002  v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
6003 
6004  if (WantTwoLevelState)
6005  {
6006 #ifdef MFEM_USE_MEMALLOC
6007  BisectedElement *aux = BEMemory.Alloc();
6008  aux->SetCoarseElem(tri);
6009 #else
6010  BisectedElement *aux = new BisectedElement(tri);
6011 #endif
6012  aux->FirstChild = tri = new Triangle(v[0], tri->GetAttribute());
6013  aux->SecondChild = NumOfElements;
6014  pce[0] = aux;
6015  }
6016  else
6017  tri->SetVertices(v[0]); // changes vert[0..2] !!!
6018  elements.Append(new Triangle(v[1], tri->GetAttribute()));
6019 
6020  // 3. edge1 and edge2 may have to be changed for the second triangle.
6021  if (v[1][0] < v_to_v.NumberOfRows() && v[1][1] < v_to_v.NumberOfRows())
6022  {
6023  bisect = v_to_v(v[1][0], v[1][1]);
6024 #ifdef MFEM_DEBUG
6025  if (bisect < 0)
6026  mfem_error("Mesh::Bisection(...) of triangle! #2");
6027 #endif
6028  if (edge1[bisect] == i)
6029  edge1[bisect] = NumOfElements;
6030  else if (edge2[bisect] == i)
6031  edge2[bisect] = NumOfElements;
6032  }
6033  NumOfElements++;
6034  }
6035  else if (t == Element::TETRAHEDRON)
6036  {
6037  int j, type, new_type, old_redges[2], new_redges[2][2], flag;
6038  Tetrahedron *tet = (Tetrahedron *) pce[0];
6039 
6040  if (tet->GetRefinementFlag() == 0)
6041  mfem_error("Mesh::Bisection : TETRAHEDRON element is not marked for "
6042  "refinement.");
6043 
6044  vert = tet->GetVertices();
6045 
6046  // 1. Get the index for the new vertex in v_new.
6047  bisect = v_to_v(vert[0], vert[1]);
6048  if (bisect == -1)
6049  {
6050  tet->ParseRefinementFlag(old_redges, type, flag);
6051  cerr << "Error in Bisection(...) of tetrahedron!" << endl
6052  << " redge[0] = " << old_redges[0]
6053  << " redge[1] = " << old_redges[1]
6054  << " type = " << type
6055  << " flag = " << flag << endl;
6056  mfem_error();
6057  }
6058 
6059  if (middle[bisect] == -1)
6060  {
6061  v_new = NumOfVertices++;
6062  for (j = 0; j < 3; j++)
6063  V(j) = 0.5 * (vertices[vert[0]](j) + vertices[vert[1]](j));
6064  vertices.Append(V);
6065 
6066  middle[bisect] = v_new;
6067  }
6068  else
6069  v_new = middle[bisect];
6070 
6071  // 2. Set the node indices for the new elements in v[2][4] so that
6072  // the edge marked for refinement is between the first two nodes.
6073  tet->ParseRefinementFlag(old_redges, type, flag);
6074 
6075  v[0][3] = v_new;
6076  v[1][3] = v_new;
6077  new_redges[0][0] = 2;
6078  new_redges[0][1] = 1;
6079  new_redges[1][0] = 2;
6080  new_redges[1][1] = 1;
6081  switch (old_redges[0])
6082  {
6083  case 2:
6084  v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
6085  if (type == Tetrahedron::TYPE_PF) new_redges[0][1] = 4;
6086  break;
6087  case 3:
6088  v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
6089  break;
6090  case 5:
6091  v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
6092  }
6093  switch (old_redges[1])
6094  {
6095  case 1:
6096  v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
6097  if (type == Tetrahedron::TYPE_PF) new_redges[1][0] = 3;
6098  break;
6099  case 4:
6100  v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
6101  break;
6102  case 5:
6103  v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
6104  }
6105 
6106  int attr = tet->GetAttribute();
6107  if (WantTwoLevelState)
6108  {
6109 #ifdef MFEM_USE_MEMALLOC
6110  BisectedElement *aux = BEMemory.Alloc();
6111  aux->SetCoarseElem(tet);
6112  tet = TetMemory.Alloc();
6113  tet->SetVertices(v[0]);
6114  tet->SetAttribute(attr);
6115 #else
6116  BisectedElement *aux = new BisectedElement(tet);
6117  tet = new Tetrahedron(v[0], attr);
6118 #endif
6119  aux->FirstChild = tet;
6120  aux->SecondChild = NumOfElements;
6121  pce[0] = aux;
6122  // 'tet' now points to the first child
6123  }
6124  else
6125  tet->SetVertices(v[0]);
6126 
6127  {
6128 #ifdef MFEM_USE_MEMALLOC
6129  Tetrahedron *tet2 = TetMemory.Alloc();
6130  tet2->SetVertices(v[1]);
6131  tet2->SetAttribute(attr);
6132  elements.Append(tet2);
6133 #else
6134  elements.Append(new Tetrahedron(v[1], attr));
6135 #endif
6136  }
6137 
6138  // 3. Set the bisection flag
6139  switch (type)
6140  {
6141  case Tetrahedron::TYPE_PU:
6142  new_type = Tetrahedron::TYPE_PF; break;
6143  case Tetrahedron::TYPE_PF:
6144  new_type = Tetrahedron::TYPE_A; break;
6145  default:
6146  new_type = Tetrahedron::TYPE_PU;
6147  }
6148 
6149  tet->CreateRefinementFlag(new_redges[0], new_type, flag+1);
6150  ((Tetrahedron *)elements[NumOfElements])->
6151  CreateRefinementFlag(new_redges[1], new_type, flag+1);
6152 
6153  NumOfElements++;
6154  }
6155  else
6156  mfem_error("Bisection for now works only for triangles & tetrahedra.");
6157 }
6158 
6159 void Mesh::Bisection(int i, const DSTable &v_to_v, int *middle)
6160 {
6161  int *vert;
6162  int v[2][3], v_new, bisect, t;
6163  Element **pce = &boundary[i];
6164 
6165  t = pce[0]->GetType();
6166  if (WantTwoLevelState)
6167  {
6168  while (1)
6169  {
6170  if (t == Element::BISECTED)
6171  pce = & ( ((BisectedElement *) pce[0])->FirstChild );
6172  else if (t == Element::QUADRISECTED)
6173  pce = & ( ((QuadrisectedElement *) pce[0])->FirstChild );
6174  else
6175  break;
6176  t = pce[0]->GetType();
6177  }
6178  }
6179 
6180  if (t == Element::TRIANGLE)
6181  {
6182  Triangle *tri = (Triangle *) pce[0];
6183 
6184  vert = tri->GetVertices();
6185 
6186  // 1. Get the index for the new vertex in v_new.
6187  bisect = v_to_v(vert[0], vert[1]);
6188 #ifdef MFEM_DEBUG
6189  if (bisect < 0)
6190  mfem_error("Mesh::Bisection(...) of boundary triangle! #1");
6191 #endif
6192  v_new = middle[bisect];
6193 #ifdef MFEM_DEBUG
6194  if (v_new == -1)
6195  mfem_error("Mesh::Bisection(...) of boundary triangle! #2");
6196 #endif
6197 
6198  // 2. Set the node indices for the new elements in v[0] and v[1] so that
6199  // the edge marked for refinement is between the first two nodes.
6200  v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
6201  v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
6202  if (WantTwoLevelState)
6203  {
6204 #ifdef MFEM_USE_MEMALLOC
6205  BisectedElement *aux = BEMemory.Alloc();
6206  aux->SetCoarseElem(tri);
6207 #else
6208  BisectedElement *aux = new BisectedElement(tri);
6209 #endif
6210  aux->FirstChild = tri = new Triangle(v[0], tri->GetAttribute());
6211  aux->SecondChild = NumOfBdrElements;
6212  pce[0] = aux;
6213  // 'tri' now points to the first child
6214  }
6215  else
6216  tri->SetVertices(v[0]);
6217  boundary.Append(new Triangle(v[1], tri->GetAttribute()));
6218 
6219  NumOfBdrElements++;
6220  }
6221  else
6222  mfem_error("Bisection of boundary elements works only for triangles!");
6223 }
6224 
6225 void Mesh::UniformRefinement(int i, const DSTable &v_to_v,
6226  int *edge1, int *edge2, int *middle)
6227 {
6228  Array<int> v;
6229  int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
6230  Vertex V;
6231 
6232  if (elements[i]->GetType() == Element::TRIANGLE)
6233  {
6234  elements[i]->GetVertices(v);
6235 
6236  // 1. Get the indeces for the new vertices in array v_new
6237  bisect[0] = v_to_v(v[0],v[1]);
6238  bisect[1] = v_to_v(v[1],v[2]);
6239  bisect[2] = v_to_v(v[0],v[2]);
6240 #ifdef MFEM_DEBUG
6241  if (bisect[0] < 0 || bisect[1] < 0 || bisect[2] < 0)
6242  mfem_error("Mesh::UniformRefinement(...): ERROR");
6243 #endif
6244 
6245  for (j = 0; j < 3; j++) // for the 3 edges fix v_new
6246  if (middle[bisect[j]] == -1)
6247  {
6248  v_new[j] = NumOfVertices++;
6249  for (int d = 0; d < spaceDim; d++)
6250  V(d) = (vertices[v[j]](d) + vertices[v[(j+1)%3]](d))/2.;
6251  vertices.Append(V);
6252 
6253  // Put the element that may need refinement (because of this
6254  // bisection) in edge1, or -1 if no more refinement is needed.
6255  if (edge1[bisect[j]] == i)
6256  edge1[bisect[j]] = edge2[bisect[j]];
6257 
6258  middle[bisect[j]] = v_new[j];
6259  }
6260  else
6261  {
6262  v_new[j] = middle[bisect[j]];
6263 
6264  // This edge will require no more refinement.
6265  edge1[bisect[j]] = -1;
6266  }
6267 
6268  // 2. Set the node indeces for the new elements in v1, v2, v3 & v4 so that
6269  // the edges marked for refinement be between the first two nodes.
6270  v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
6271  v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
6272  v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
6273  v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
6274 
6275  elements.Append(new Triangle(v1, elements[i]->GetAttribute()));
6276  elements.Append(new Triangle(v2, elements[i]->GetAttribute()));
6277  elements.Append(new Triangle(v3, elements[i]->GetAttribute()));
6278  if (WantTwoLevelState)
6279  {
6280  QuadrisectedElement *aux = new QuadrisectedElement(elements[i]);
6281  aux->FirstChild = new Triangle(v4, elements[i]->GetAttribute());
6282  aux->Child2 = NumOfElements;
6283  aux->Child3 = NumOfElements+1;
6284  aux->Child4 = NumOfElements+2;
6285  elements[i] = aux;
6286  }
6287  else
6288  {
6289  elements[i]->SetVertices(v4);
6290  }
6291 
6292  NumOfElements += 3;
6293  }
6294  else
6295  mfem_error("Uniform refinement for now works only for triangles.");
6296 }
6297 
6298 void Mesh::SetState(int s)
6299 {
6300  if (ncmesh)
6301  {
6302  if (State != Mesh::NORMAL && s == Mesh::NORMAL)
6303  {
6304  delete nc_coarse_level;
6305  nc_coarse_level = NULL;
6306  ncmesh->ClearCoarseLevel();
6307  }
6308  else if ((State == Mesh::TWO_LEVEL_COARSE && s == Mesh::TWO_LEVEL_FINE) ||
6309  (State == Mesh::TWO_LEVEL_FINE && s == Mesh::TWO_LEVEL_COARSE))
6310  {
6311  this->Swap(*nc_coarse_level, false);
6312  State = s;
6313  }
6314  else if (State != s)
6315  {
6316  mfem_error("Oops! Mesh::SetState");
6317  }
6318  }
6319  else
6320  {
6321  if (State != Mesh::NORMAL && s == Mesh::NORMAL)
6322  {
6323  // two level state --->> normal state
6324  int i, t;
6325 
6326  for (i = 0; i < f_NumOfElements; )
6327  {
6328  t = elements[i]->GetType();
6329  if (t == Element::BISECTED ||
6330  t == Element::QUADRISECTED ||
6331  t == Element::OCTASECTED)
6332  {
6333  RefinedElement *aux = (RefinedElement *) elements[i];
6334  elements[i] = aux->FirstChild;
6335  FreeElement(aux->CoarseElem);
6336  FreeElement(aux);
6337  }
6338  else
6339  i++;
6340  }
6341 
6342  for (i = 0; i < f_NumOfBdrElements; )
6343  {
6344  t = boundary[i]->GetType();
6345  if (t == Element::BISECTED ||
6346  t == Element::QUADRISECTED ||
6347  t == Element::OCTASECTED)
6348  {
6349  RefinedElement *aux = (RefinedElement *) boundary[i];
6350  boundary[i] = aux->FirstChild;
6351  FreeElement(aux->CoarseElem);
6352  FreeElement(aux);
6353  }
6354  else
6355  i++;
6356  }
6357 
6358  if (el_to_edge != NULL)
6359  {
6360  delete c_el_to_edge;
6361  el_to_edge = f_el_to_edge;
6362  if (Dim == 2)
6363  {
6364  if (State == Mesh::TWO_LEVEL_COARSE)
6365  mfem::Swap(be_to_edge, fc_be_to_edge);
6366  fc_be_to_edge.DeleteAll();
6367  }
6368  if (Dim == 3)
6369  {
6370  delete c_bel_to_edge;
6371  bel_to_edge = f_bel_to_edge;
6372  }
6373  }
6374  if (el_to_face != NULL)
6375  {
6376  delete c_el_to_face;
6377  el_to_face = f_el_to_face;
6378  if (State == Mesh::TWO_LEVEL_COARSE)
6379  mfem::Swap(faces_info, fc_faces_info);
6380  fc_faces_info.DeleteAll();
6381  }
6382 
6383  NumOfVertices = f_NumOfVertices;
6384  NumOfEdges = f_NumOfEdges;
6385  if (Dim == 3)
6386  NumOfFaces = f_NumOfFaces;
6387  NumOfElements = f_NumOfElements;
6388  NumOfBdrElements = f_NumOfBdrElements;
6390  State = s;
6391  }
6392  else if (State == Mesh::TWO_LEVEL_COARSE && s == Mesh::TWO_LEVEL_FINE)
6393  {
6394  if (el_to_edge != NULL)
6395  {
6396  el_to_edge = f_el_to_edge;
6397  if (Dim == 2)
6398  mfem::Swap(be_to_edge, fc_be_to_edge);
6399  if (Dim == 3)
6400  bel_to_edge = f_bel_to_edge;
6401  }
6402  if (el_to_face != NULL)
6403  {
6404  el_to_face = f_el_to_face;
6405  mfem::Swap(faces_info, fc_faces_info);
6406  }
6407  NumOfVertices = f_NumOfVertices;
6408  NumOfEdges = f_NumOfEdges;
6409  if (Dim == 3)
6410  NumOfFaces = f_NumOfFaces;
6411  NumOfElements = f_NumOfElements;
6412  NumOfBdrElements = f_NumOfBdrElements;
6414  State = s;
6415  }
6416  else if (State == Mesh::TWO_LEVEL_FINE && s == Mesh::TWO_LEVEL_COARSE)
6417  {
6418  if (el_to_edge != NULL)
6419  {
6420  el_to_edge = c_el_to_edge;
6421  if (Dim == 2)
6422  mfem::Swap(be_to_edge, fc_be_to_edge);
6423  if (Dim == 3)
6424  bel_to_edge = c_bel_to_edge;
6425  }
6426  if (el_to_face != NULL)
6427  {
6428  el_to_face = c_el_to_face;
6429  mfem::Swap(faces_info, fc_faces_info);
6430  }
6431  NumOfVertices = c_NumOfVertices;
6432  NumOfEdges = c_NumOfEdges;
6433  if (Dim == 3)
6434  NumOfFaces = c_NumOfFaces;
6435  NumOfElements = c_NumOfElements;
6436  NumOfBdrElements = c_NumOfBdrElements;
6438  State = s;
6439  }
6440  else if (State != s)
6441  mfem_error("Oops! Mesh::SetState");
6442  }
6443 }
6444 
6446 {
6447  int t = elements[i]->GetType();
6448 
6449  if (Dim == 1)
6450  {
6451  if (t == Element::BISECTED)
6452  return 2;
6453  }
6454  else if (Dim == 2)
6455  {
6456  if (t == Element::QUADRISECTED)
6457  return 4;
6458  else if (t == Element::BISECTED)
6459  {
6460  // assuming that the elements are either BisectedElements or
6461  // regular elements
6462  int n = 1;
6463  BisectedElement *aux = (BisectedElement *) elements[i];
6464  do
6465  {
6466  n += GetNumFineElems(aux->SecondChild);
6467  if (aux->FirstChild->GetType() != Element::BISECTED)
6468  break;
6469  aux = (BisectedElement *) (aux->FirstChild);
6470  }
6471  while (1);
6472  return n;
6473  }
6474  }
6475  else if (Dim == 3)
6476  {
6477  // assuming that the element is a BisectedElement,
6478  // OctasectedElement (with children that are regular elements) or
6479  // regular element
6480  if (t == Element::BISECTED)
6481  {
6482  int n = 1;
6483  BisectedElement *aux = (BisectedElement *) elements[i];
6484  do
6485  {
6486  n += GetNumFineElems (aux->SecondChild);
6487  if (aux->FirstChild->GetType() != Element::BISECTED)
6488  break;
6489  aux = (BisectedElement *) (aux->FirstChild);
6490  }
6491  while (1);
6492  return n;
6493  }
6494  else if (t == Element::OCTASECTED)
6495  return 8;
6496  return 1; // regular element (i.e. it is not refined)
6497  }
6498 
6499  return 1; // the element is not refined
6500 }
6501 
6503 {
6504  if (E->GetType() == Element::BISECTED)
6505  {
6506  int L, R, n, s, lb, rb;
6507 
6508  L = GetBisectionHierarchy(((BisectedElement *)E)->FirstChild);
6509  n = ((BisectedElement *)E)->SecondChild;
6510  R = GetBisectionHierarchy(elements[n]);
6511  n = 1; s = 1;
6512  lb = rb = 1;
6513  do
6514  {
6515  int nlb, nrb;
6516  nlb = nrb = 0;
6517  while (lb > 0)
6518  {
6519  n |= ((L & 1) << s);
6520  s++;
6521  nlb += (L & 1);
6522  L = (L >> 1);
6523  lb--;
6524  }
6525  while (rb > 0)
6526  {
6527  n |= ((R & 1) << s);
6528  s++;
6529  nrb += (R & 1);
6530  R = (R >> 1);
6531  rb--;
6532  }
6533  lb = 2 * nlb; rb = 2 * nrb;
6534  }
6535  while (lb > 0 || rb > 0);
6536  return n;
6537  }
6538  return 0;
6539 }
6540 
6542 {
6543  int t;
6544 
6545  if (Dim == 1)
6546  {
6547  if (elements[i]->GetType() == Element::BISECTED)
6548  return 1; // refinement type for bisected SEGMENT
6549  }
6550  else if (Dim == 2)
6551  {
6552  t = elements[i]->GetType();
6553  if (t == Element::QUADRISECTED)
6554  {
6555  t = ((QuadrisectedElement *)elements[i])->CoarseElem->GetType();
6556  if (t == Element::QUADRILATERAL)
6557  return 1; // refinement type for quadrisected QUADRILATERAL
6558  else
6559  return 2; // refinement type for quadrisected TRIANGLE
6560  }
6561  else if (t == Element::BISECTED)
6562  {
6563  int type;
6564  type = GetBisectionHierarchy(elements[i]);
6565  if (type == 0)
6566  mfem_error("Mesh::GetRefinementType(...)");
6567  return type+2;
6568  }
6569  }
6570  else if (Dim == 3)
6571  {
6572  int redges[2], type, flag;
6573  Element *E = elements[i];
6574  Tetrahedron *tet;
6575 
6576  t = E->GetType();
6577  if (t != Element::BISECTED)
6578  {
6579  if (t == Element::OCTASECTED)
6580  return 1; // refinement type for octasected CUBE
6581  else
6582  return 0;
6583  }
6584  // Bisected TETRAHEDRON
6585  tet = (Tetrahedron *) (((BisectedElement *) E)->CoarseElem);
6586  tet->ParseRefinementFlag(redges, type, flag);
6587  if (type == Tetrahedron::TYPE_A && redges[0] == 2)
6588  type = 5;
6589  else if (type == Tetrahedron::TYPE_M && redges[0] == 2)
6590  type = 6;
6591  type++;
6592  type |= ( GetBisectionHierarchy(E) << 3 );
6593  if (type < 8) type = 0;
6594 
6595  return type;
6596  }
6597 
6598  return 0; // no refinement
6599 }
6600 
6601 int Mesh::GetFineElem(int i, int j)
6602 {
6603  int t = elements[i]->GetType();
6604 
6605  if (Dim == 1)
6606  {
6607  if (t == Element::BISECTED)
6608  {
6609  switch (j)
6610  {
6611  case 0: return i;
6612  default: return ((BisectedElement *)elements[i])->SecondChild;
6613  }
6614  }
6615  }
6616  else if (Dim == 2)
6617  {
6618  if (t == Element::QUADRISECTED)
6619  {
6620  QuadrisectedElement *aux = (QuadrisectedElement *) elements[i];
6621  if (aux->CoarseElem->GetType() == Element::QUADRILATERAL)
6622  switch (j)
6623  {
6624  case 0: return i;
6625  case 1: return aux->Child2;
6626  case 2: return aux->Child3;
6627  case 3: return aux->Child4;
6628  default: mfem_error("Mesh::GetFineElem #1");
6629  }
6630  else // quadrisected TRIANGLE
6631  switch (j)
6632  {
6633  case 0: return aux->Child2;
6634  case 1: return aux->Child3;
6635  case 2: return aux->Child4;
6636  case 3: return i;
6637  default: mfem_error("Mesh::GetFineElem #2");
6638  }
6639  }
6640  else if (t == Element::BISECTED)
6641  {
6642  int n = 0;
6643  BisectedElement *aux = (BisectedElement *) elements[i];
6644  do
6645  {
6646  int k = GetFineElem(aux->SecondChild, j-n);
6647  if (k >= 0)
6648  return k;
6649  n -= k; // (-k) is the number of the leaves in this SecondChild
6650  // n is the number of the leaves in
6651  // the SecondChild-ren so far
6652  if (aux->FirstChild->GetType() != Element::BISECTED)
6653  break;
6654  aux = (BisectedElement *) (aux->FirstChild);
6655  }
6656  while (1);
6657  if (j > n) // i.e. if (j >= n+1)
6658  return -(n+1);
6659  return i; // j == n, i.e. j is the index of the last leaf
6660  }
6661  }
6662  else if (Dim == 3)
6663  {
6664  if (t == Element::BISECTED)
6665  {
6666  int n = 0;
6667  BisectedElement *aux = (BisectedElement *) elements[i];
6668  do
6669  {
6670  int k = GetFineElem(aux->SecondChild, j-n);
6671  if (k >= 0)
6672  return k;
6673  n -= k; // (-k) is the number of the leaves in this SecondChild
6674  // n is the number of the leaves in
6675  // the SecondChild-ren so far
6676  if (aux->FirstChild->GetType() != Element::BISECTED)
6677  break;
6678  aux = (BisectedElement *) (aux->FirstChild);
6679  }
6680  while (1);
6681  if (j > n) // i.e. if (j >= n+1)
6682  return -(n+1);
6683  return i; // j == n, i.e. j is the index of the last leaf
6684  }
6685  else if (t == Element::OCTASECTED)
6686  {
6687  if (j == 0) return i;
6688  return ((OctasectedElement *) elements[i])->Child[j-1];
6689  }
6690  }
6691 
6692  if (j > 0)
6693  return -1;
6694 
6695  return i; // no refinement
6696 }
6697 
6698 void Mesh::BisectTriTrans(DenseMatrix &pointmat, Triangle *tri, int child)
6699 {
6700  double np[2];
6701 
6702  if (child == 0) // left triangle
6703  {
6704  // Set the new coordinates of the vertices
6705  np[0] = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
6706  np[1] = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
6707  pointmat(0,1) = pointmat(0,0); pointmat(1,1) = pointmat(1,0);
6708  pointmat(0,0) = pointmat(0,2); pointmat(1,0) = pointmat(1,2);
6709  pointmat(0,2) = np[0]; pointmat(1,2) = np[1];
6710  }
6711  else // right triangle
6712  {
6713  // Set the new coordinates of the vertices
6714  np[0] = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
6715  np[1] = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
6716  pointmat(0,0) = pointmat(0,1); pointmat(1,0) = pointmat(1,1);
6717  pointmat(0,1) = pointmat(0,2); pointmat(1,1) = pointmat(1,2);
6718  pointmat(0,2) = np[0]; pointmat(1,2) = np[1];
6719  }
6720 }
6721 
6722 void Mesh::BisectTetTrans(DenseMatrix &pointmat, Tetrahedron *tet, int child)
6723 {
6724  int i, j, redges[2], type, flag, ind[4];
6725  double t[4];
6726 
6727  tet->ParseRefinementFlag(redges, type, flag);
6728 
6729  if (child == 0) // left tetrahedron
6730  {
6731  // Set the new coordinates of the vertices
6732  pointmat(0,1) = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
6733  pointmat(1,1) = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
6734  pointmat(2,1) = 0.5 * ( pointmat(2,0) + pointmat(2,1) );
6735  // Permute the vertices according to the type & redges
6736  switch (redges[0])
6737  {
6738  case 2: ind[0] = 0; ind[1] = 3; ind[2] = 1; ind[3] = 2; break;
6739  case 3: ind[0] = 1; ind[1] = 3; ind[2] = 2; ind[3] = 0; break;
6740  case 5:
6741  default: ind[0] = 2; ind[1] = 3; ind[2] = 0; ind[3] = 1;
6742  }
6743  }
6744  else // right tetrahedron
6745  {
6746  // Set the new coordinates of the vertices
6747  pointmat(0,0) = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
6748  pointmat(1,0) = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
6749  pointmat(2,0) = 0.5 * ( pointmat(2,0) + pointmat(2,1) );
6750  // Permute the vertices according to the type & redges
6751  switch (redges[1])
6752  {
6753  case 1: ind[0] = 3; ind[1] = 1; ind[2] = 0; ind[3] = 2; break;
6754  case 4: ind[0] = 3; ind[1] = 0; ind[2] = 2; ind[3] = 1; break;
6755  case 5:
6756  default: ind[0] = 3; ind[1] = 2; ind[2] = 1; ind[3] = 0;
6757  }
6758  }
6759  // Do the permutation
6760  for (i = 0; i < 3; i++)
6761  {
6762  for (j = 0; j < 4; j++)
6763  t[j] = pointmat(i,j);
6764  for (j = 0; j < 4; j++)
6765  pointmat(i,ind[j]) = t[j];
6766  }
6767 }
6768 
6769 int Mesh::GetFineElemPath(int i, int j)
6770 {
6771  // if (Dim == 3)
6772  {
6773  if (elements[i]->GetType() == Element::BISECTED)
6774  {
6775  int n = 0, l = 0;
6776  BisectedElement *aux = (BisectedElement *) elements[i];
6777  do
6778  {
6779  int k = GetFineElemPath(aux->SecondChild, j-n);
6780  if (k >= 0)
6781  return ((k << 1)+1) << l;
6782  n -= k; // (-k) is the number of the leaves in this SecondChild
6783  // n is the number of the leaves in
6784  // the SecondChild-ren so far
6785  l++;
6786  if (aux->FirstChild->GetType() != Element::BISECTED)
6787  break;
6788  aux = (BisectedElement *) (aux->FirstChild);
6789  }
6790  while (1);
6791  if (j > n) // i.e. if (j >= n+1)
6792  return -(n+1);
6793  return 0; // j == n, i.e. j is the index of the last leaf
6794  }
6795  if (j > 0)
6796  return -1;
6797  }
6798 
6799  return 0;
6800 }
6801 
6803 {
6804  int t;
6805 
6806  if (Dim == 1)
6807  {
6808  DenseMatrix &pm = Transformation.GetPointMat();
6809  Transformation.SetFE(&SegmentFE);
6810  Transformation.Attribute = 0;
6811  Transformation.ElementNo = 0;
6812  pm.SetSize(1, 2);
6813  if (elements[i]->GetType() == Element::BISECTED)
6814  {
6815  switch (j)
6816  {
6817  case 0: pm(0,0) = 0.0; pm(0,1) = 0.5; break;
6818  default: pm(0,0) = 0.5; pm(0,1) = 1.0; break;
6819  }
6820  }
6821  else
6822  {
6823  pm(0,0) = 0.0; pm(0,1) = 1.0;
6824  }
6825  }
6826  else if (Dim == 2)
6827  {
6828  DenseMatrix &pm = Transformation.GetPointMat();
6829  Transformation.Attribute = 0;
6830  Transformation.ElementNo = 0;
6831  t = elements[i]->GetType();
6832  if (t == Element::QUADRISECTED)
6833  {
6834  t = ((QuadrisectedElement *)elements[i])->CoarseElem->GetType();
6835  if (t == Element::QUADRILATERAL)
6836  {
6837  // quadrisected QUADRILATERAL
6838  Transformation.SetFE(&QuadrilateralFE);
6839  pm.SetSize(2, 4);
6840  switch (j)
6841  {
6842  case 0:
6843  pm(0,0) = 0.0; pm(1,0) = 0.0; // x; y;
6844  pm(0,1) = 0.5; pm(1,1) = 0.0;
6845  pm(0,2) = 0.5; pm(1,2) = 0.5;
6846  pm(0,3) = 0.0; pm(1,3) = 0.5;
6847  break;
6848  case 1:
6849  pm(0,0) = 0.5; pm(1,0) = 0.0;
6850  pm(0,1) = 1.0; pm(1,1) = 0.0;
6851  pm(0,2) = 1.0; pm(1,2) = 0.5;
6852  pm(0,3) = 0.5; pm(1,3) = 0.5;
6853  break;
6854  case 2:
6855  pm(0,0) = 0.5; pm(1,0) = 0.5;
6856  pm(0,1) = 1.0; pm(1,1) = 0.5;
6857  pm(0,2) = 1.0; pm(1,2) = 1.0;
6858  pm(0,3) = 0.5; pm(1,3) = 1.0;
6859  break;
6860  case 3:
6861  pm(0,0) = 0.0; pm(1,0) = 0.5;
6862  pm(0,1) = 0.5; pm(1,1) = 0.5;
6863  pm(0,2) = 0.5; pm(1,2) = 1.0;
6864  pm(0,3) = 0.0; pm(1,3) = 1.0;
6865  break;
6866  default:
6867  mfem_error("Mesh::GetFineElemTrans(...) 1");
6868  }
6869  }
6870  else
6871  {
6872  // quadrisected TRIANGLE
6873  Transformation.SetFE(&TriangleFE);
6874  pm.SetSize(2, 3);
6875  switch (j)
6876  {
6877  case 0:
6878  pm(0,0) = 0.0; pm(0,1) = 0.5; pm(0,2) = 0.0; // x
6879  pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 0.5; // y
6880  break;
6881  case 1:
6882  pm(0,0) = 0.5; pm(0,1) = 1.0; pm(0,2) = 0.5;
6883  pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 0.5;
6884  break;
6885  case 2:
6886  pm(0,0) = 0.0; pm(0,1) = 0.5; pm(0,2) = 0.0;
6887  pm(1,0) = 0.5; pm(1,1) = 0.5; pm(1,2) = 1.0;
6888  break;
6889  case 3:
6890  pm(0,0) = 0.5; pm(0,1) = 0.0; pm(0,2) = 0.5;
6891  pm(1,0) = 0.5; pm(1,1) = 0.5; pm(1,2) = 0.0;
6892  break;
6893  default:
6894  mfem_error("Mesh::GetFineElemTrans(...) 2");
6895  }
6896  }
6897  }
6898  else if (t == Element::BISECTED)
6899  {
6900  // bisected TRIANGLE
6901  Transformation.SetFE(&TriangleFE);
6902  pm.SetSize(2, 3);
6903 
6904  int path;
6905  Element *E;
6906 
6907  // pm is initialzed with the coordinates of the vertices of the
6908  // reference triangle
6909  pm(0,0) = 0.0; pm(0,1) = 1.0; pm(0,2) = 0.0;
6910  pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 1.0;
6911 
6912  path = GetFineElemPath(i, j);
6913 
6914  E = elements[i];
6915  while (E->GetType() == Element::BISECTED)
6916  {
6917  BisectedElement *aux = (BisectedElement *) E;
6918 
6919  BisectTriTrans(pm, (Triangle *) aux->CoarseElem, path & 1);
6920  E = (path & 1) ? elements[aux->SecondChild] : aux->FirstChild;
6921  path = path >> 1;
6922  }
6923  }
6924  else
6925  {
6926  // identity transformation
6927  Transformation.SetFE(&TriangleFE);
6928  pm.SetSize(2, 3);
6929  pm(0,0) = 0.0; pm(0,1) = 1.0; pm(0,2) = 0.0;
6930  pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 1.0;
6931  }
6932  return &Transformation;
6933  }
6934  else if (Dim == 3)
6935  {
6936  if (elements[i]->GetType() == Element::OCTASECTED)
6937  {
6938  int jj;
6939  double dx, dy, dz;
6940  DenseMatrix &pm = Transformation.GetPointMat();
6941  Transformation.SetFE(&HexahedronFE);
6942  Transformation.Attribute = 0;
6943  Transformation.ElementNo = 0;
6944