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