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