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