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