phynexis-deps/mfem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
29ae8fb6ac9fce927fdfd01ef796f9df38f76e97
mfem/mesh/mesh.cpp
T
Tzanio Kolev 29ae8fb6ac Version 1.2
2011-04-08 15:27:24 -07:00

6614 lines
180 KiB
C++
Raw Blame History

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.googlecode.com.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
// Implementation of data type mesh
#include <iostream>
#include <fstream>
#include <math.h>
#include <string.h>
#include <time.h>
#include "mesh_headers.hpp"
#include "../fem/fem.hpp"
#include "../general/sort_pairs.hpp"
void Mesh::GetElementJacobian(int i, DenseMatrix &J)
{
int geom = GetElementBaseGeometry(i);
ElementTransformation *eltransf = GetElementTransformation(i);
eltransf->SetIntPoint(&Geometries.GetCenter(geom));
Geometries.JacToPerfJac(geom, eltransf->Jacobian(), J);
}
double Mesh::GetElementSize(int i, int type)
{
DenseMatrix J(Dim);
GetElementJacobian(i, J);
if (type)
{
Vector sv;
J.SingularValues(sv);
if (type == 1)
return sv(sv.Size()-1); // h_min
return sv(0); // h_max
}
return pow(fabs(J.Det()), 1./Dim);
}
double Mesh::GetElementSize(int i, const Vector &dir)
{
DenseMatrix J(Dim);
Vector d_hat(Dim);
GetElementJacobian(i, J);
J.MultTranspose(dir, d_hat);
return sqrt((d_hat * d_hat) / (dir * dir));
}
double Mesh::GetElementVolume(int i)
{
ElementTransformation *et = GetElementTransformation(i);
const IntegrationRule &ir = IntRules.Get(GetElementBaseGeometry(i),
et->OrderJ());
double volume = 0.0;
for (int j = 0; j < ir.GetNPoints(); j++)
{
const IntegrationPoint &ip = ir.IntPoint(j);
et->SetIntPoint(&ip);
volume += ip.weight * et->Weight();
}
return volume;
}
void Mesh::PrintCharacteristics(Vector *Vh, Vector *Vk)
{
int i, dim;
DenseMatrix J;
double h_min, h_max, kappa_min, kappa_max, h, kappa;
#ifdef MFEM_USE_LAPACK
Vector sv;
#else
DenseMatrix Jinv;
#endif
cout << "Mesh Characteristics:" << flush;
dim = Dimension();
J.SetSize(dim);
#ifndef MFEM_USE_LAPACK
Jinv.SetSize(dim);
#endif
if (Vh) Vh->SetSize(NumOfElements);
if (Vk) Vk->SetSize(NumOfElements);
for (i = 0; i < NumOfElements; i++)
{
GetElementJacobian(i, J);
h = pow(fabs(J.Det()), 1.0/double(dim));
#ifdef MFEM_USE_LAPACK
// J's condition number in spectral norm
J.SingularValues(sv);
kappa = sv(0) / sv(dim-1);
#else
// J's condition number in Frobenius norm
CalcInverse(J, Jinv);
kappa = J.FNorm() * Jinv.FNorm();
#endif
if (Vh) (*Vh)(i) = h;
if (Vk) (*Vk)(i) = kappa;
if (i == 0)
{
h_min = h_max = h;
kappa_min = kappa_max = kappa;
}
else
{
if (h < h_min) h_min = h;
if (h > h_max) h_max = h;
if (kappa < kappa_min) kappa_min = kappa;
if (kappa > kappa_max) kappa_max = kappa;
}
}
if (dim == 2)
cout << endl
<< "Number of vertices : " << GetNV() << endl
<< "Number of edges : " << GetNEdges() << endl
<< "Number of elements : " << GetNE() << endl
<< "Number of bdr elem : " << GetNBE() << endl
<< "Euler Number : " << EulerNumber2D() << endl
<< "h_min : " << h_min << endl
<< "h_max : " << h_max << endl
<< "kappa_min : " << kappa_min << endl
<< "kappa_max : " << kappa_max << endl
<< endl;
else
cout << endl
<< "Number of vertices : " << GetNV() << endl
<< "Number of edges : " << GetNEdges() << endl
<< "Number of faces : " << GetNFaces() << endl
<< "Number of elements : " << GetNE() << endl
<< "Number of bdr elem : " << GetNBE() << endl
<< "Euler Number : " << EulerNumber() << endl
<< "h_min : " << h_min << endl
<< "h_max : " << h_max << endl
<< "kappa_min : " << kappa_min << endl
<< "kappa_max : " << kappa_max << endl
<< endl;
}
FiniteElement *Mesh::GetTransformationFEforElementType(int ElemType)
{
switch (ElemType)
{
case Element::POINT : return &PointFE;
case Element::SEGMENT : return &SegmentFE;
case Element::TRIANGLE : return &TriangleFE;
case Element::QUADRILATERAL : return &QuadrilateralFE;
case Element::TETRAHEDRON : return &TetrahedronFE;
case Element::HEXAHEDRON : return &HexahedronFE;
}
mfem_error("Mesh::GetTransformationFEforElement - unknown ElementType");
return &TriangleFE;
}
void Mesh::GetElementTransformation(int i, IsoparametricTransformation *ElTr)
{
ElTr->Attribute = GetAttribute(i);
ElTr->ElementNo = i;
if (Nodes == NULL)
{
GetPointMatrix(i, ElTr->GetPointMat());
ElTr->SetFE(GetTransformationFEforElementType(GetElementType(i)));
}
else
{
DenseMatrix &pm = ElTr->GetPointMat();
Array<int> vdofs;
Nodes->FESpace()->GetElementVDofs(i, vdofs);
int n = vdofs.Size()/Dim;
pm.SetSize(Dim, n);
for (int k = 0; k < Dim; k++)
for (int j = 0; j < n; j++)
pm(k,j) = (*Nodes)(vdofs[n*k+j]);
ElTr->SetFE(Nodes->FESpace()->GetFE(i));
}
}
void Mesh::GetElementTransformation(int i, const Vector &nodes,
IsoparametricTransformation *ElTr)
{
ElTr->Attribute = GetAttribute(i);
ElTr->ElementNo = i;
DenseMatrix &pm = ElTr->GetPointMat();
if (Nodes == NULL)
{
int nv = elements[i]->GetNVertices();
const int *v = elements[i]->GetVertices();
int n = vertices.Size();
pm.SetSize(Dim, nv);
for (int k = 0; k < Dim; k++)
for (int j = 0; j < nv; j++)
pm(k, j) = nodes(k*n+v[j]);
ElTr->SetFE(GetTransformationFEforElementType(GetElementType(i)));
}
else
{
Array<int> vdofs;
Nodes->FESpace()->GetElementVDofs(i, vdofs);
int n = vdofs.Size()/Dim;
pm.SetSize(Dim, n);
for (int k = 0; k < Dim; k++)
for (int j = 0; j < n; j++)
pm(k,j) = nodes(vdofs[n*k+j]);
ElTr->SetFE(Nodes->FESpace()->GetFE(i));
}
}
ElementTransformation *Mesh::GetElementTransformation(int i)
{
GetElementTransformation(i, &Transformation);
return &Transformation;
}
ElementTransformation *Mesh::GetBdrElementTransformation(int i)
{
FaceTransformation.Attribute = GetBdrAttribute(i);
FaceTransformation.ElementNo = i; // boundary element number
if (Nodes == NULL)
{
GetBdrPointMatrix(i, FaceTransformation.GetPointMat());
FaceTransformation.SetFE(
GetTransformationFEforElementType(GetBdrElementType(i)));
}
else
{
DenseMatrix &pm = FaceTransformation.GetPointMat();
Array<int> vdofs;
Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
int n = vdofs.Size()/Dim;
pm.SetSize(Dim, n);
for (int k = 0; k < Dim; k++)
for (int j = 0; j < n; j++)
pm(k,j) = (*Nodes)(vdofs[n*k+j]);
FaceTransformation.SetFE(Nodes->FESpace()->GetBE(i));
}
return &FaceTransformation;
}
void Mesh::GetLocalSegToTriTransformation(
IsoparametricTransformation &Transf, int i)
{
static const int tri_faces[3][2] = {{1, 0}, {2, 1}, {0, 2}};
static const int seg_inv_orient[2][2] = {{0, 1}, {1, 0}};
int j;
const int *tv, *so;
const IntegrationRule *TriVert;
DenseMatrix &locpm = Transf.GetPointMat();
Transf.SetFE(&SegmentFE);
tv = tri_faces[i/64]; // (i/64) is the local face no. in the triangle
so = seg_inv_orient[i%64]; // (i%64) is the orientation of the segment
TriVert = Geometries.GetVertices(Geometry::TRIANGLE);
locpm.SetSize(2, 2);
for (j = 0; j < 2; j++)
{
locpm(0, so[j]) = TriVert->IntPoint(tv[j]).x;
locpm(1, so[j]) = TriVert->IntPoint(tv[j]).y;
}
}
void Mesh::GetLocalSegToQuadTransformation(
IsoparametricTransformation &Transf, int i)
{
static const int quad_faces[4][2] = {{1, 0}, {2, 1}, {3, 2}, {0, 3}};
static const int seg_inv_orient[2][2] = {{0, 1}, {1, 0}};
int j;
const int *qv, *so;
const IntegrationRule *QuadVert;
DenseMatrix &locpm = Transf.GetPointMat();
Transf.SetFE(&SegmentFE);
qv = quad_faces[i/64]; // (i/64) is the local face no. in the quad
so = seg_inv_orient[i%64]; // (i%64) is the orientation of the segment
QuadVert = Geometries.GetVertices(Geometry::SQUARE);
locpm.SetSize(2, 2);
for (j = 0; j < 2; j++)
{
locpm(0, so[j]) = QuadVert->IntPoint(qv[j]).x;
locpm(1, so[j]) = QuadVert->IntPoint(qv[j]).y;
}
}
void Mesh::GetLocalTriToTetTransformation(
IsoparametricTransformation &Transf, int i)
{
static const int tet_faces[4][3] = {{1, 2, 3}, {0, 3, 2},
{0, 1, 3}, {0, 2, 1}};
static const int tri_inv_orient[6][3] = {{0, 1, 2}, {1, 0, 2},
{1, 2, 0}, {2, 1, 0},
{2, 0, 1}, {0, 2, 1}};
int j;
const int *tv, *to;
const IntegrationRule *TetVert;
DenseMatrix &locpm = Transf.GetPointMat();
Transf.SetFE(&TriangleFE);
tv = tet_faces[i/64]; // (i/64) is the local face no. in the tet
// (i%64) is the orientation of the tetrahedron face
// w.r.t. the face element
to = tri_inv_orient[i%64];
TetVert = Geometries.GetVertices(Geometry::TETRAHEDRON);
locpm.SetSize(3, 3);
for (j = 0; j < 3; j++)
{
locpm(0, to[j]) = TetVert->IntPoint(tv[j]).x;
locpm(1, to[j]) = TetVert->IntPoint(tv[j]).y;
locpm(2, to[j]) = TetVert->IntPoint(tv[j]).z;
}
}
void Mesh::GetLocalQuadToHexTransformation(
IsoparametricTransformation &Transf, int i)
{
static const int hex_faces[6][4] = {{3, 2, 1, 0}, {0, 1, 5, 4},
{1, 2, 6, 5}, {2, 3, 7, 6},
{3, 0, 4, 7}, {4, 5, 6, 7}};
// must be 'quad_inv_or' ... fix me
static const int quad_orient[8][4] = {{0, 1, 2, 3}, {0, 3, 2, 1},
{1, 2, 3, 0}, {1, 0, 3, 2},
{2, 3, 0, 1}, {2, 1, 0, 3},
{3, 0, 1, 2}, {3, 2, 1, 0}};
int j;
const int *hv, *qo;
const IntegrationRule *HexVert;
DenseMatrix &locpm = Transf.GetPointMat();
Transf.SetFE(&QuadrilateralFE);
hv = hex_faces[i/64]; // (i/64) is the local face no. in the hex
qo = quad_orient[i%64]; // (i%64) is the orientation of the quad
HexVert = Geometries.GetVertices(Geometry::CUBE);
locpm.SetSize(3, 4);
for (j = 0; j < 4; j++)
{
locpm(0, qo[j]) = HexVert->IntPoint(hv[j]).x;
locpm(1, qo[j]) = HexVert->IntPoint(hv[j]).y;
locpm(2, qo[j]) = HexVert->IntPoint(hv[j]).z;
}
}
FaceElementTransformations *Mesh::GetFaceElementTransformations(int FaceNo)
{
int i, j, *v, nv;
FiniteElement *FaceFE;
// setup the transformation for the first element
FaceElemTr.Elem1No = faces_info[FaceNo].Elem1No;
GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
FaceElemTr.Elem1 = &Transformation;
// setup the transformation for the second element
// return NULL in the Elem2 field if there's no second element, i.e.
// the face is on the "boundary"
if ( (FaceElemTr.Elem2No = faces_info[FaceNo].Elem2No) >= 0 )
{
GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
FaceElemTr.Elem2 = &Transformation2;
}
else
FaceElemTr.Elem2 = NULL;
FaceElemTr.FaceGeom = faces[FaceNo]->GetGeometryType();
FaceFE = GetTransformationFEforElementType(faces[FaceNo]->GetType());
// setup the face transformation
FaceTransformation.Attribute = faces[FaceNo]->GetAttribute();
FaceTransformation.ElementNo = FaceNo;
v = faces[FaceNo]->GetVertices();
nv = faces[FaceNo]->GetNVertices();
DenseMatrix &pm = FaceTransformation.GetPointMat();
pm.SetSize(Dim, nv);
for (i = 0; i < Dim; i++)
for (j = 0; j < nv; j++)
pm(i, j) = vertices[v[j]](i);
FaceTransformation.SetFE(FaceFE);
FaceElemTr.Face = &FaceTransformation;
// setup Loc1 & Loc2
switch (faces[FaceNo]->GetType())
{
case Element::SEGMENT:
if (GetElementType(faces_info[FaceNo].Elem1No) == Element::TRIANGLE)
GetLocalSegToTriTransformation(FaceElemTr.Loc1.Transf,
faces_info[FaceNo].Elem1Inf);
else // assume the element is a quad
GetLocalSegToQuadTransformation(FaceElemTr.Loc1.Transf,
faces_info[FaceNo].Elem1Inf);
if (FaceElemTr.Elem2No >= 0)
if (GetElementType(faces_info[FaceNo].Elem2No)
== Element::TRIANGLE)
GetLocalSegToTriTransformation(FaceElemTr.Loc2.Transf,
faces_info[FaceNo].Elem2Inf);
else // assume the element is a quad
GetLocalSegToQuadTransformation(FaceElemTr.Loc2.Transf,
faces_info[FaceNo].Elem2Inf);
break;
case Element::TRIANGLE:
// --------- assumes the face is a triangle -- face of a tetrahedron
GetLocalTriToTetTransformation(FaceElemTr.Loc1.Transf,
faces_info[FaceNo].Elem1Inf);
if (FaceElemTr.Elem2No >= 0)
GetLocalTriToTetTransformation(FaceElemTr.Loc2.Transf,
faces_info[FaceNo].Elem2Inf);
break;
case Element::QUADRILATERAL:
// --------- assumes the face is a quad -- face of a hexahedron
GetLocalQuadToHexTransformation(FaceElemTr.Loc1.Transf,
faces_info[FaceNo].Elem1Inf);
if (FaceElemTr.Elem2No >= 0)
GetLocalQuadToHexTransformation(FaceElemTr.Loc2.Transf,
faces_info[FaceNo].Elem2Inf);
break;
}
return &FaceElemTr;
}
FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
{
FaceElementTransformations *tr;
int fn;
if (Dim == 3)
fn = be_to_face[BdrElemNo];
else
fn = be_to_edge[BdrElemNo];
if (faces_info[fn].Elem2No >= 0)
return NULL;
tr = GetFaceElementTransformations(fn);
tr->Face->Attribute = boundary[BdrElemNo]->GetAttribute();
return tr;
}
void Mesh::GetFaceElements(int Face, int *Elem1, int *Elem2)
{
*Elem1 = faces_info[Face].Elem1No;
*Elem2 = faces_info[Face].Elem2No;
}
void Mesh::GetFaceInfos(int Face, int *Inf1, int *Inf2)
{
*Inf1 = faces_info[Face].Elem1Inf;
*Inf2 = faces_info[Face].Elem2Inf;
}
void Mesh::Init()
{
NumOfVertices = NumOfElements = NumOfBdrElements = NumOfEdges = -1;
WantTwoLevelState = 0;
State = Mesh::NORMAL;
Nodes = NULL;
own_nodes = 1;
}
void Mesh::InitTables()
{
el_to_edge = el_to_face = el_to_el =
bel_to_edge = face_edge = edge_vertex = NULL;
}
void Mesh::DeleteTables()
{
if (el_to_edge != NULL)
delete el_to_edge;
if (el_to_face != NULL)
delete el_to_face;
if (el_to_el != NULL)
delete el_to_el;
if (Dim == 3 && bel_to_edge != NULL)
delete bel_to_edge;
if (face_edge != NULL)
delete face_edge;
if (edge_vertex != NULL)
delete edge_vertex;
InitTables();
}
void Mesh::DeleteCoarseTables()
{
delete el_to_el;
delete face_edge;
delete edge_vertex;
el_to_el = face_edge = edge_vertex = NULL;
}
void Mesh::SetAttributes()
{
int i, j, nattr;
Array<int> attribs;
attribs.SetSize(GetNBE());
for (i = 0; i < attribs.Size(); i++)
attribs[i] = GetBdrAttribute(i);
attribs.Sort();
if (attribs.Size() > 0)
nattr = 1;
else
nattr = 0;
for (i = 1; i < attribs.Size(); i++)
if (attribs[i] != attribs[i-1])
nattr++;
bdr_attributes.SetSize(nattr);
if (nattr > 0)
{
bdr_attributes[0] = attribs[0];
for (i = j = 1; i < attribs.Size(); i++)
if (attribs[i] != attribs[i-1])
bdr_attributes[j++] = attribs[i];
if (attribs[0] <= 0)
cout << "Mesh::SetAttributes(): "
"Non-positive attributes on the boundary!"
<< endl;
}
attribs.SetSize(GetNE());
for (i = 0; i < attribs.Size(); i++)
attribs[i] = GetAttribute(i);
attribs.Sort();
if (attribs.Size() > 0)
nattr = 1;
else
nattr = 0;
for (i = 1; i < attribs.Size(); i++)
if (attribs[i] != attribs[i-1])
nattr++;
attributes.SetSize(nattr);
if (nattr > 0)
{
attributes[0] = attribs[0];
for (i = j = 1; i < attribs.Size(); i++)
if (attribs[i] != attribs[i-1])
attributes[j++] = attribs[i];
if (attribs[0] <= 0)
cout << "Mesh::SetAttributes(): "
"Non-positive attributes in the domain!"
<< endl;
}
}
Mesh::Mesh(int _Dim, int NVert, int NElem, int NBdrElem)
{
Dim = _Dim;
Init();
InitTables();
NumOfVertices = 0;
vertices.SetSize(NVert); // just allocate space for vertices
NumOfElements = 0;
elements.SetSize(NElem); // just allocate space for Element *
NumOfBdrElements = 0;
boundary.SetSize(NBdrElem); // just allocate space for Element *
}
void Mesh::AddVertex(double *x)
{
double *y = vertices[NumOfVertices]();
for (int i = 0; i < Dim; i++)
y[i] = x[i];
NumOfVertices++;
}
void Mesh::AddTri(int *vi, int attr)
{
elements[NumOfElements++] = new Triangle(vi, attr);
}
void Mesh::AddTriangle(int *vi, int attr)
{
elements[NumOfElements++] = new Triangle(vi, attr);
}
void Mesh::AddQuad(int *vi, int attr)
{
elements[NumOfElements++] = new Quadrilateral(vi, attr);
}
void Mesh::AddTet(int *vi, int attr)
{
#ifdef MFEM_USE_MEMALLOC
Tetrahedron *tet;
tet = TetMemory.Alloc();
tet->SetVertices(vi);
tet->SetAttribute(attr);
elements[NumOfElements++] = tet;
#else
elements[NumOfElements++] = new Tetrahedron(vi, attr);
#endif
}
void Mesh::AddHex(int *vi, int attr)
{
elements[NumOfElements++] = new Hexahedron(vi, attr);
}
void Mesh::AddBdrSegment(int *vi, int attr)
{
boundary[NumOfBdrElements++] = new Segment(vi, attr);
}
void Mesh::AddBdrTriangle(int *vi, int attr)
{
boundary[NumOfBdrElements++] = new Triangle(vi, attr);
}
void Mesh::AddBdrQuad(int *vi, int attr)
{
boundary[NumOfBdrElements++] = new Quadrilateral(vi, attr);
}
void Mesh::GenerateBoundaryElements()
{
int i, j;
// GenerateFaces();
for (i = 0; i < boundary.Size(); i++)
FreeElement(boundary[i]);
// count the 'NumOfBdrElements'
NumOfBdrElements = 0;
for (i = 0; i < faces_info.Size(); i++)
if (faces_info[i].Elem2No == -1)
NumOfBdrElements++;
boundary.SetSize(NumOfBdrElements);
if (Dim == 3)
{
be_to_face.SetSize(NumOfBdrElements);
delete bel_to_edge;
bel_to_edge = NULL;
}
if (Dim == 2)
be_to_edge.SetSize(NumOfBdrElements);
for (j = i = 0; i < faces_info.Size(); i++)
if (faces_info[i].Elem2No == -1)
{
boundary[j] = faces[i]->Duplicate();
if (Dim == 3)
be_to_face[j] = i;
if (Dim == 2)
be_to_edge[j] = i;
j++;
}
// In 3D, 'bel_to_edge' is destroyed but it's not updated.
}
typedef struct {
int edge;
double length;
} edge_length;
// Used by qsort to sort edges in increasing (according their length) order.
static int edge_compare(const void *ii, const void *jj)
{
edge_length *i = (edge_length *)ii, *j = (edge_length *)jj;
if (i->length > j->length) return (1);
if (i->length < j->length) return (-1);
return (0);
}
void Mesh::FinalizeTriMesh(int generate_edges, int refine)
{
CheckElementOrientation();
if (refine)
MarkTriMeshForRefinement();
if (generate_edges)
{
el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
GenerateFaces();
CheckBdrElementOrientation();
}
else
NumOfEdges = 0;
NumOfFaces = 0;
SetAttributes();
meshgen = 1;
}
void Mesh::FinalizeQuadMesh(int generate_edges, int refine)
{
CheckElementOrientation();
if (generate_edges)
{
el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
GenerateFaces();
CheckBdrElementOrientation();
}
else
NumOfEdges = 0;
NumOfFaces = 0;
SetAttributes();
meshgen = 2;
}
void Mesh::MarkForRefinement()
{
if (meshgen & 1)
{
if (Dim == 2)
MarkTriMeshForRefinement();
else if (Dim == 3)
MarkTetMeshForRefinement();
}
}
void Mesh::MarkTriMeshForRefinement()
{
// Mark the longest triangle edge by rotating the indeces so that
// vertex 0 - vertex 1 to be the longest element's edge.
DenseMatrix pmat;
for (int i = 0; i < NumOfElements; i++)
if (elements[i]->GetType() == Element::TRIANGLE)
{
GetPointMatrix(i, pmat);
elements[i]->MarkEdge(pmat);
}
}
void Mesh::MarkTetMeshForRefinement()
{
// Mark the longest tetrahedral edge by rotating the indices so that
// vertex 0 - vertex 1 is the longest edge in the element.
DSTable v_to_v(NumOfVertices);
GetVertexToVertexTable(v_to_v);
NumOfEdges = v_to_v.NumberOfEntries();
edge_length *length = new edge_length[NumOfEdges];
for (int i = 0; i < NumOfVertices; i++)
{
for (DSTable::RowIterator it(v_to_v, i); !it; ++it)
{
int j = it.Index();
length[j].length = GetLength(i, it.Column());
length[j].edge = j;
}
}
// sort in increasing order
qsort(length, NumOfEdges, sizeof(edge_length), edge_compare);
int *order = new int [NumOfEdges];
for (int i = 0; i < NumOfEdges; i++)
order[length[i].edge] = i;
for (int i = 0; i < NumOfElements; i++)
if (elements[i]->GetType() == Element::TETRAHEDRON)
elements[i]->MarkEdge(v_to_v, order);
for (int i = 0; i < NumOfBdrElements; i++)
if (boundary[i]->GetType() == Element::TRIANGLE)
boundary[i]->MarkEdge(v_to_v, order);
delete [] order;
delete [] length;
}
void Mesh::FinalizeTetMesh(int generate_edges, int refine)
{
CheckElementOrientation();
if (refine)
{
MarkTetMeshForRefinement();
}
GetElementToFaceTable();
GenerateFaces();
if (NumOfBdrElements == 0)
GenerateBoundaryElements();
if (generate_edges == 1)
{
el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
}
else
{
el_to_edge = NULL; // Not really necessary -- InitTables was called
bel_to_edge = NULL;
NumOfEdges = 0;
}
SetAttributes();
meshgen = 1;
}
void Mesh::FinalizeHexMesh(int generate_edges, int refine)
{
CheckElementOrientation();
GetElementToFaceTable();
GenerateFaces();
if (NumOfBdrElements == 0)
GenerateBoundaryElements();
CheckBdrElementOrientation();
if (generate_edges)
{
el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
}
else
NumOfEdges = 0;
SetAttributes();
meshgen = 2;
}
Mesh::Mesh(int nx, int ny, Element::Type type, int generate_edges,
double sx, double sy)
{
int i, j, k;
Dim = 2;
Init();
InitTables();
// Creates quadrilateral mesh
if (type == Element::QUADRILATERAL)
{
meshgen = 2;
NumOfVertices = (nx+1) * (ny+1);
NumOfElements = nx * ny;
NumOfBdrElements = 2 * nx + 2 * ny;
vertices.SetSize(NumOfVertices);
elements.SetSize(NumOfElements);
boundary.SetSize(NumOfBdrElements);
double cx, cy;
int ind[4];
// Sets vertices and the corresponding coordinates
k = 0;
for (j = 0; j < ny+1; j++)
{
cy = ((double) j / ny) * sy;
for (i = 0; i < nx+1; i++)
{
cx = ((double) i / nx) * sx;
vertices[k](0) = cx;
vertices[k](1) = cy;
k++;
}
}
// Sets elements and the corresponding indices of vertices
k = 0;
for (j = 0; j < ny; j++)
{
for (i = 0; i < nx; i++)
{
ind[0] = i + j*(nx+1);
ind[1] = i + 1 +j*(nx+1);
ind[2] = i + 1 + (j+1)*(nx+1);
ind[3] = i + (j+1)*(nx+1);
elements[k] = new Quadrilateral(ind);
k++;
}
}
// Sets boundary elements and the corresponding indices of vertices
int m = (nx+1)*ny;
for (i = 0; i < nx; i++)
{
boundary[i] = new Segment(i, i+1, 1);
boundary[nx+i] = new Segment(m+i, m+i+1, 3);
}
m = nx+1;
for (j = 0; j < ny; j++)
{
boundary[2*nx+j] = new Segment(j*m, (j+1)*m, 4);
boundary[2*nx+ny+j] = new Segment(j*m+nx, (j+1)*m+nx, 2);
}
}
// Creates triangular mesh
if (type == Element::TRIANGLE)
{
meshgen = 1;
NumOfVertices = (nx+1) * (ny+1);
NumOfElements = 2 * nx * ny;
NumOfBdrElements = 2 * nx + 2 * ny;
vertices.SetSize(NumOfVertices);
elements.SetSize(NumOfElements);
boundary.SetSize(NumOfBdrElements);
double cx, cy;
int ind[3];
// Sets vertices and the corresponding coordinates
k = 0;
for (j = 0; j < ny+1; j++)
{
cy = ((double) j / ny) * sy;
for (i = 0; i < nx+1; i++)
{
cx = ((double) i / nx) * sx;
vertices[k](0) = cx;
vertices[k](1) = cy;
k++;
}
}
// Sets the elements and the corresponding indices of vertices
k = 0;
for (j = 0; j < ny; j++)
{
for (i = 0; i < nx; i++)
{
ind[0] = i + j*(nx+1);
ind[1] = i + 1 + (j+1)*(nx+1);
ind[2] = i + (j+1)*(nx+1);
elements[k] = new Triangle(ind);
k++;
ind[1] = i + 1 + j*(nx+1);
ind[2] = i + 1 + (j+1)*(nx+1);
elements[k] = new Triangle(ind);
k++;
}
}
// Sets boundary elements and the corresponding indices of vertices
int m = (nx+1)*ny;
for (i = 0; i < nx; i++)
{
boundary[i] = new Segment(i, i+1, 1);
boundary[nx+i] = new Segment(m+i, m+i+1, 3);
}
m = nx+1;
for (j = 0; j < ny; j++)
{
boundary[2*nx+j] = new Segment(j*m, (j+1)*m, 4);
boundary[2*nx+ny+j] = new Segment(j*m+nx, (j+1)*m+nx, 2);
}
MarkTriMeshForRefinement();
}
CheckElementOrientation();
if (generate_edges == 1)
{
el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
GenerateFaces();
CheckBdrElementOrientation();
}
else
NumOfEdges = 0;
NumOfFaces = 0;
attributes.Append(1);
bdr_attributes.Append(1); bdr_attributes.Append(2);
bdr_attributes.Append(3); bdr_attributes.Append(4);
}
Mesh::Mesh(int n)
{
int j, ind[1];
Dim = 1;
Init();
InitTables();
meshgen = 1;
NumOfVertices = n + 1;
NumOfElements = n;
NumOfBdrElements = 2;
vertices.SetSize(NumOfVertices);
elements.SetSize(NumOfElements);
boundary.SetSize(NumOfBdrElements);
// Sets vertices and the corresponding coordinates
for (j = 0; j < n+1; j++)
vertices[j](0) = (double) j / n;
// Sets elements and the corresponding indices of vertices
for (j = 0; j < n; j++)
elements[j] = new Segment(j, j+1, 1);
// Sets the boundary elements
ind[0] = 0;
boundary[0] = new Point(ind, 1);
ind[0] = n;
boundary[1] = new Point(ind, 2);
NumOfEdges = 0;
NumOfFaces = 0;
}
Mesh::Mesh(istream &input, int generate_edges, int refine)
{
Init();
InitTables();
Load(input, generate_edges, refine);
}
// mesh types:
// 0 - MFEM mesh v1.0
// 1 - 1D mesh ("linemesh"?)
// 2 - 2D Netgen mesh ("areamesh2") / MFEM curvil. ext. ("curved_areamesh2")
// 3 - 3D Netgen mesh ("NETGEN" / "NETGEN_Neutral_Format")
// 4 - TrueGrid mesh (2D/3D how?)
// 5 - VTK mesh
static int GetMeshTypeFromStream(istream &meshin)
{
string buff;
int type = -1;
streampos start_pos = meshin.tellg();
meshin >> ws;
getline(meshin, buff);
if (buff == "MFEM mesh v1.0")
type = 0;
if (buff == "linemesh")
type = 1;
else if (buff == "areamesh2" || buff == "curved_areamesh2")
type = 2;
else if (buff == "NETGEN" || buff == "NETGEN_Neutral_Format")
type = 3;
else if (buff == "TrueGrid")
type = 4;
else if (buff == "# vtk DataFile Version 3.0")
type = 5;
meshin.seekg(start_pos);
return type;
}
Element *NewElement(int geom)
{
switch (geom)
{
case Geometry::SEGMENT: return (new Segment);
case Geometry::TRIANGLE: return (new Triangle);
case Geometry::SQUARE: return (new Quadrilateral);
case Geometry::CUBE: return (new Hexahedron);
case Geometry::TETRAHEDRON:
#ifdef MFEM_USE_MEMALLOC
return TetMemory.Alloc();
#else
return (new Tetrahedron);
#endif
}
return NULL;
}
// see Tetrahedron::edges
static const int vtk_quadratic_tet[10] =
{ 0, 1, 2, 3, 4, 7, 5, 6, 8, 9 };
// see Hexahedron::edges & Mesh::GenerateFaces
static const int vtk_quadratic_hex[27] =
{ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
24, 22, 21, 23, 20, 25, 26 };
void Mesh::Load(istream &input, int generate_edges, int refine)
{
int i, j, ints[32], n, attr, curved = 0, read_gf = 1;
const int buflen = 1024;
char buf[buflen];
#ifdef MFEM_DEBUG
if (!input)
mfem_error("Input file stream not opened : Mesh::Load");
#endif
if (NumOfVertices != -1)
{
// Delete the elements.
for (i = 0; i < NumOfElements; i++)
// delete elements[i];
FreeElement(elements[i]);
elements.DeleteAll();
// Delete the vertices.
vertices.DeleteAll();
// Delete the boundary elements.
for (i = 0; i < NumOfBdrElements; i++)
// delete boundary[i];
FreeElement(boundary[i]);
boundary.DeleteAll();
// Delete interior faces (if generated)
for (i = 0; i < faces.Size(); i++)
FreeElement(faces[i]);
faces.DeleteAll();
faces_info.DeleteAll();
// Delete the edges (if generated).
DeleteTables();
be_to_edge.DeleteAll();
be_to_face.DeleteAll();
}
InitTables();
if (own_nodes) delete Nodes;
Nodes = NULL;
int mesh_type = GetMeshTypeFromStream(input);
if (mesh_type < 0)
{
mfem_error("Mesh::Load : Unknown input mesh format!");
return;
}
if (mesh_type == 0)
{
// Read MFEM mesh v1.0 format
string ident;
int attr, geom, nv, *v;
// read the mesh format line and lines begining with '#' (comments)
do
{
input >> ws;
getline(input, ident);
input >> ws;
}
while (input.peek() == '#');
input >> ident; // 'dimension'
input >> Dim;
input >> ident; // 'elements'
input >> NumOfElements;
elements.SetSize(NumOfElements);
for (j = 0; j < NumOfElements; j++)
{
input >> attr >> geom;
elements[j] = NewElement(geom);
elements[j]->SetAttribute(attr);
nv = elements[j]->GetNVertices();
v = elements[j]->GetVertices();
for (i = 0; i < nv; i++)
input >> v[i];
}
input >> ident; // 'boundary'
input >> NumOfBdrElements;
boundary.SetSize(NumOfBdrElements);
for (j = 0; j < NumOfBdrElements; j++)
{
input >> attr >> geom;
boundary[j] = NewElement(geom);
boundary[j]->SetAttribute(attr);
nv = boundary[j]->GetNVertices();
v = boundary[j]->GetVertices();
for (i = 0; i < nv; i++)
input >> v[i];
}
input >> ident; // 'vertices'
input >> NumOfVertices;
vertices.SetSize(NumOfVertices);
input >> ws >> ident;
if (ident != "nodes")
{
// read the vertices
int vdim = atoi(ident.c_str());
for (j = 0; j < NumOfVertices; j++)
for (i = 0; i < vdim; i++)
input >> vertices[j](i);
}
else
{
// prepare to read the nodes
input >> ws;
curved = 1;
}
}
else if (mesh_type == 1)
{
int j,p1,p2,a;
Dim = 1;
// Read 1st (info) line
input >> buf;
input >> NumOfVertices;
vertices.SetSize(NumOfVertices);
// Sets vertices and the corresponding coordinates
for (j = 0; j < NumOfVertices; j++)
input >> vertices[j](0);
input >> NumOfElements;
elements.SetSize(NumOfElements);
// Sets elements and the corresponding indices of vertices
for (j = 0; j < NumOfElements; j++)
{
input >> a >> p1 >> p2;
elements[j] = new Segment(p1-1, p2-1, a);
}
int ind[1];
input >> NumOfBdrElements;
boundary.SetSize(NumOfBdrElements);
for (j = 0; j < NumOfBdrElements; j++)
{
input >> a >> ind[0];
ind[0]--;
boundary[j] = new Point(ind,a);
}
}
else if (mesh_type == 2)
{
// Read planar mesh in Netgen format.
Dim = 2;
// Read the type of the mesh.
input >> buf;
if (!strcmp("curved_areamesh2", buf))
curved = 1;
// Read the boundary elements.
input >> NumOfBdrElements;
boundary.SetSize(NumOfBdrElements);
for (i = 0; i < NumOfBdrElements; i++)
{
input >> attr
>> ints[0] >> ints[1];
ints[0]--; ints[1]--;
boundary[i] = new Segment(ints, attr);
}
// Read the elements.
input >> NumOfElements;
elements.SetSize(NumOfElements);
for (i = 0; i < NumOfElements; i++)
{
input >> attr >> n;
for (j = 0; j < n; j++)
{
input >> ints[j];
ints[j]--;
}
switch (n)
{
case 2:
elements[i] = new Segment(ints, attr);
break;
case 3:
elements[i] = new Triangle(ints, attr);
break;
case 4:
elements[i] = new Quadrilateral(ints, attr);
break;
}
}
if (!curved)
{
// Read the vertices.
input >> NumOfVertices;
vertices.SetSize(NumOfVertices);
for (i = 0; i < NumOfVertices; i++)
for (j = 0; j < Dim; j++)
input >> vertices[i](j);
}
else
{
input >> NumOfVertices;
vertices.SetSize(NumOfVertices);
input >> ws;
}
}
else if (mesh_type == 3)
{
// Read a netgen format mesh of tetrahedra.
Dim = 3;
// Read the type of the mesh.
input >> buf;
// Read the vertices
input >> NumOfVertices;
vertices.SetSize(NumOfVertices);
for (i = 0; i < NumOfVertices; i++)
for (j = 0; j < Dim; j++)
input >> vertices[i](j);
// Read the elements
input >> NumOfElements;
elements.SetSize(NumOfElements);
for (i = 0; i < NumOfElements; i++)
{
input >> attr;
for (j = 0; j < 4; j++)
{
input >> ints[j];
ints[j]--;
}
#ifdef MFEM_USE_MEMALLOC
Tetrahedron *tet;
tet = TetMemory.Alloc();
tet->SetVertices(ints);
tet->SetAttribute(attr);
elements[i] = tet;
#else
elements[i] = new Tetrahedron(ints, attr);
#endif
}
// Read the boundary information.
input >> NumOfBdrElements;
boundary.SetSize(NumOfBdrElements);
for (i = 0; i < NumOfBdrElements; i++)
{
input >> attr;
for (j = 0; j < 3; j++)
{
input >> ints[j];
ints[j]--;
}
boundary[i] = new Triangle(ints, attr);
}
}
else if (mesh_type == 4)
{
// Reading TrueGrid mesh.
// Read the type of the mesh.
input >> buf;
// TODO: find the actual dimension
Dim = 3;
if (Dim == 2)
{
int vari;
double varf;
input >> vari >> NumOfVertices >> vari >> vari >> NumOfElements;
input.getline(buf, buflen);
input.getline(buf, buflen);
input >> vari;
input.getline(buf, buflen);
input.getline(buf, buflen);
input.getline(buf, buflen);
// Read the vertices.
vertices.SetSize(NumOfVertices);
for (i = 0; i < NumOfVertices; i++)
{
input >> vari >> varf >> vertices[i](0) >> vertices[i](1);
input.getline(buf, buflen);
}
// Read the elements.
elements.SetSize(NumOfElements);
for (i = 0; i < NumOfElements; i++)
{
input >> vari >> attr;
for (j = 0; j < 4; j++)
{
input >> ints[j];
ints[j]--;
}
input.getline(buf, buflen);
input.getline(buf, buflen);
elements[i] = new Quadrilateral(ints, attr);
}
}
else if (Dim == 3)
{
int vari;
double varf;
input >> vari >> NumOfVertices >> NumOfElements;
input.getline(buf, buflen);
input.getline(buf, buflen);
input >> vari >> vari >> NumOfBdrElements;
input.getline(buf, buflen);
input.getline(buf, buflen);
input.getline(buf, buflen);
// Read the vertices.
vertices.SetSize(NumOfVertices);
for (i = 0; i < NumOfVertices; i++)
{
input >> vari >> varf >> vertices[i](0) >> vertices[i](1)
>> vertices[i](2);
input.getline(buf, buflen);
}
// Read the elements.
elements.SetSize(NumOfElements);
for (i = 0; i < NumOfElements; i++)
{
input >> vari >> attr;
for (j = 0; j < 8; j++)
{
input >> ints[j];
ints[j]--;
}
input.getline(buf, buflen);
elements[i] = new Hexahedron(ints, attr);
}
// Read the boundary elements.
boundary.SetSize(NumOfBdrElements);
for (i = 0; i < NumOfBdrElements; i++)
{
input >> attr;
for (j = 0; j < 4; j++)
{
input >> ints[j];
ints[j]--;
}
input.getline(buf, buflen);
boundary[i] = new Quadrilateral(ints, attr);
}
}
}
else if (mesh_type == 5)
{
// Reading VTK mesh
string buff;
getline(input, buff); // "# vtk DataFile Version 3.0"
getline(input, buff); // comment line
getline(input, buff);
if (buff != "ASCII")
{
mfem_error("Mesh::Load : VTK mesh is not in ASCII format!");
return;
}
getline(input, buff);
if (buff != "DATASET UNSTRUCTURED_GRID")
{
mfem_error("Mesh::Load : VTK mesh is not UNSTRUCTURED_GRID!");
return;
}
// Read the points
int np = 0;
Vector points;
input >> buff;
if (buff == "POINTS")
{
input >> np >> ws;
points.SetSize(3*np);
getline(input, buff); // "double"
for (i = 0; i < points.Size(); i++)
input >> points(i);
}
// Read the cells
NumOfElements = n = 0;
Array<int> cells_data;
input >> ws >> buff;
if (buff == "CELLS")
{
input >> NumOfElements >> n >> ws;
cells_data.SetSize(n);
for (i = 0; i < n; i++)
input >> cells_data[i];
}
// Read the cell types
Dim = 0;
int order = 1;
input >> ws >> buff;
if (buff == "CELL_TYPES")
{
input >> NumOfElements;
elements.SetSize(NumOfElements);
for (j = i = 0; i < NumOfElements; i++)
{
int ct;
input >> ct;
switch (ct)
{
case 5: // triangle
Dim = 2;
elements[i] = new Triangle(&cells_data[j+1]);
break;
case 9: // quadrilateral
Dim = 2;
elements[i] = new Quadrilateral(&cells_data[j+1]);
break;
case 10: // tetrahedron
Dim = 3;
#ifdef MFEM_USE_MEMALLOC
elements[i] = TetMemory.Alloc();
elements[i]->SetVertices(&cells_data[j+1]);
#else
elements[i] = new Tetrahedron(&cells_data[j+1]);
#endif
break;
case 12: // hexahedron
Dim = 3;
elements[i] = new Hexahedron(&cells_data[j+1]);
break;
case 22: // quadratic triangle
Dim = 2;
order = 2;
elements[i] = new Triangle(&cells_data[j+1]);
break;
case 28: // biquadratic quadrilateral
Dim = 2;
order = 2;
elements[i] = new Quadrilateral(&cells_data[j+1]);
break;
case 24: // quadratic tetrahedron
Dim = 3;
order = 2;
#ifdef MFEM_USE_MEMALLOC
elements[i] = TetMemory.Alloc();
elements[i]->SetVertices(&cells_data[j+1]);
#else
elements[i] = new Tetrahedron(&cells_data[j+1]);
#endif
break;
case 29: // triquadratic hexahedron
Dim = 3;
order = 2;
elements[i] = new Hexahedron(&cells_data[j+1]);
break;
default:
cerr << "Mesh::Load : VTK mesh : cell type " << ct
<< " is not supported!" << endl;
mfem_error();
return;
}
j += cells_data[j] + 1;
}
}
// Read attributes
streampos sp = input.tellg();
input >> ws >> buff;
if (buff == "CELL_DATA")
{
input >> n >> ws;
getline(input, buff);
if (buff == "SCALARS material int")
{
getline(input, buff); // "LOOKUP_TABLE default"
for (i = 0; i < NumOfElements; i++)
{
input >> attr;
elements[i]->SetAttribute(attr);
}
}
else
input.seekg(sp);
}
else
input.seekg(sp);
if (order == 1)
{
cells_data.DeleteAll();
NumOfVertices = np;
vertices.SetSize(np);
for (i = 0; i < np; i++)
{
vertices[i](0) = points(3*i+0);
vertices[i](1) = points(3*i+1);
vertices[i](2) = points(3*i+2);
}
points.Destroy();
// No boundary is defined in a VTK mesh
NumOfBdrElements = 0;
}
else if (order == 2)
{
curved = 1;
// generate new enumeration for the vertices
Array<int> pts_dof(np);
pts_dof = -1;
for (n = i = 0; i < NumOfElements; i++)
{
int *v = elements[i]->GetVertices();
int nv = elements[i]->GetNVertices();
for (j = 0; j < nv; j++)
if (pts_dof[v[j]] == -1)
pts_dof[v[j]] = n++;
}
// keep the original ordering of the vertices
for (n = i = 0; i < np; i++)
if (pts_dof[i] != -1)
pts_dof[i] = n++;
// update the element vertices
for (i = 0; i < NumOfElements; i++)
{
int *v = elements[i]->GetVertices();
int nv = elements[i]->GetNVertices();
for (j = 0; j < nv; j++)
v[j] = pts_dof[v[j]];
}
// Define the 'vertices' from the 'points' through the 'pts_dof' map
NumOfVertices = n;
vertices.SetSize(n);
for (i = 0; i < np; i++)
{
if ((j = pts_dof[i]) != -1)
{
vertices[j](0) = points(3*i+0);
vertices[j](1) = points(3*i+1);
vertices[j](2) = points(3*i+2);
}
}
// No boundary is defined in a VTK mesh
NumOfBdrElements = 0;
// Generate faces and edges so that we can define quadratic
// FE space on the mesh
// Generate faces
if (Dim > 2)
{
GetElementToFaceTable();
GenerateFaces();
}
else
NumOfFaces = 0;
// Generate edges
el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
if (Dim == 2)
GenerateFaces(); // 'Faces' in 2D refers to the edges
// Define quadratic FE space
FiniteElementCollection *fec = new QuadraticFECollection;
FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim);
Nodes = new GridFunction(fes);
Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
own_nodes = 1;
// Map vtk points to edge/face/element dofs
Array<int> dofs;
for (n = i = 0; i < NumOfElements; i++)
{
fes->GetElementDofs(i, dofs);
const int *vtk_mfem;
switch (elements[i]->GetGeometryType())
{
case Geometry::TRIANGLE:
case Geometry::SQUARE:
vtk_mfem = vtk_quadratic_hex; break; // identity map
case Geometry::TETRAHEDRON:
vtk_mfem = vtk_quadratic_tet; break;
case Geometry::CUBE:
vtk_mfem = vtk_quadratic_hex; break;
}
for (n++, j = 0; j < dofs.Size(); j++, n++)
{
if (pts_dof[cells_data[n]] == -1)
{
pts_dof[cells_data[n]] = dofs[vtk_mfem[j]];
}
else
{
if (pts_dof[cells_data[n]] != dofs[vtk_mfem[j]])
mfem_error("Mesh::Load : VTK mesh : "
"inconsistent quadratic mesh!");
}
}
}
// Define the 'Nodes' from the 'points' through the 'pts_dof' map
for (i = 0; i < np; i++)
{
dofs.SetSize(1);
if ((dofs[0] = pts_dof[i]) != -1)
{
fes->DofsToVDofs(dofs);
for (j = 0; j < dofs.Size(); j++)
(*Nodes)(dofs[j]) = points(3*i+j);
}
}
read_gf = 0;
}
}
// at this point the following should be defined:
// 1) Dim
// 2) NumOfElements, elements
// 3) NumOfBdrElements, boundary
// 4) NumOfVertices, with allocated space in vertices
// 5) curved
// 5a) if curved == 0, vertices must be defined
// 5b) if curved != 0 and read_gf != 0,
// 'input' must point to a GridFunction
// 5c) if curved != 0 and read_gf == 0,
// vertices and Nodes must be defined
// set the mesh type ('meshgen')
meshgen = 0;
for (i = 0; i < NumOfElements; i++)
{
switch (elements[i]->GetType())
{
case Element::SEGMENT:
case Element::TRIANGLE:
case Element::TETRAHEDRON:
meshgen |= 1; break;
case Element::QUADRILATERAL:
case Element::HEXAHEDRON:
meshgen |= 2;
}
}
if (!curved)
{
// check and fix element orientation
CheckElementOrientation();
if (refine)
MarkForRefinement();
}
// generate the faces
if (Dim > 2)
{
GetElementToFaceTable();
GenerateFaces();
if (NumOfBdrElements == 0)
GenerateBoundaryElements();
// check and fix boundary element orientation
if ( !(curved && (meshgen & 1)) )
CheckBdrElementOrientation();
}
else
NumOfFaces = 0;
// generate edges if requested
if (Dim > 1 && generate_edges == 1)
{
el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
if (Dim == 2)
{
GenerateFaces(); // 'Faces' in 2D refers to the edges
if (NumOfBdrElements == 0)
GenerateBoundaryElements();
// check and fix boundary element orientation
if ( !(curved && (meshgen & 1)) )
CheckBdrElementOrientation();
}
c_el_to_edge = NULL;
}
else
NumOfEdges = 0;
// generate the arrays 'attributes' and ' bdr_attributes'
SetAttributes();
if (curved)
{
if (read_gf)
{
Nodes = new GridFunction(this, input);
own_nodes = 1;
for (i = 0; i < Nodes->VectorDim(); i++)
{
Vector vert_val;
Nodes->GetNodalValues(vert_val, i+1);
for (j = 0; j < NumOfVertices; j++)
vertices[j](i) = vert_val(j);
}
}
// Check orientation and mark edges; only for triangles / tets
if (meshgen & 1)
{
FiniteElementSpace *fes = Nodes->FESpace();
const FiniteElementCollection *fec = fes->FEColl();
int num_edge_dofs = fec->DofForGeometry(Geometry::SEGMENT);
DSTable *old_v_to_v = NULL;
if (num_edge_dofs)
{
old_v_to_v = new DSTable(NumOfVertices);
GetVertexToVertexTable(*old_v_to_v);
}
// assuming all faces have the same geometry
int num_face_dofs =
(Dim < 3) ? 0 : fec->DofForGeometry(GetFaceBaseGeometry(0));
// assuming all elements have the same geometry
int num_elem_dofs = fec->DofForGeometry(GetElementBaseGeometry(0));
// check orientation and mark for refinement using just vertices
// (i.e. higher order curvature is not used)
CheckElementOrientation();
if (refine)
MarkForRefinement(); // changes topology!
// reorder the Nodes
Vector onodes = *Nodes;
Array<int> old_dofs, new_dofs;
int offset;
#ifdef MFEM_DEBUG
int redges = 0;
#endif
// vertex dofs do not need to be moved
offset = NumOfVertices * fec->DofForGeometry(Geometry::POINT);
// edge dofs:
// edge enumeration may be different but edge orientation is
// the same
if (num_edge_dofs > 0)
{
DSTable new_v_to_v(NumOfVertices);
GetVertexToVertexTable(new_v_to_v);
for (i = 0; i < NumOfVertices; i++)
{
for (DSTable::RowIterator it(new_v_to_v, i); !it; ++it)
{
int old_i = (*old_v_to_v)(i, it.Column());
int new_i = it.Index();
#ifdef MFEM_DEBUG
if (old_i != new_i)
redges++;
#endif
old_dofs.SetSize(num_edge_dofs);
new_dofs.SetSize(num_edge_dofs);
for (j = 0; j < num_edge_dofs; j++)
{
old_dofs[j] = offset + old_i * num_edge_dofs + j;
new_dofs[j] = offset + new_i * num_edge_dofs + j;
}
fes->DofsToVDofs(old_dofs);
fes->DofsToVDofs(new_dofs);
for (j = 0; j < old_dofs.Size(); j++)
(*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
}
}
offset += NumOfEdges * num_edge_dofs;
delete old_v_to_v;
}
#ifdef MFEM_DEBUG
cout << "Mesh::Load : redges = " << redges << endl;
#endif
// face dofs:
// both enumeration and orientation of the faces
// may be different
if (num_face_dofs > 0)
{
// generate the old face-vertex table
Table old_face_vertex;
old_face_vertex.MakeI(NumOfFaces);
for (i = 0; i < NumOfFaces; i++)
old_face_vertex.AddColumnsInRow(i, faces[i]->GetNVertices());
old_face_vertex.MakeJ();
for (i = 0; i < NumOfFaces; i++)
old_face_vertex.AddConnections(i, faces[i]->GetVertices(),
faces[i]->GetNVertices());
old_face_vertex.ShiftUpI();
// update 'el_to_face', 'be_to_face', 'faces', and 'faces_info'
STable3D *faces_tbl = GetElementToFaceTable(1);
GenerateFaces();
// loop over the old face numbers
for (i = 0; i < NumOfFaces; i++)
{
int *old_v = old_face_vertex.GetRow(i), *new_v;
int new_i, new_or, *dof_ord;
switch (old_face_vertex.RowSize(i))
{
case 3:
new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
new_v = faces[new_i]->GetVertices();
new_or = GetTriOrientation(old_v, new_v);
dof_ord = fec->DofOrderForOrientation(Geometry::TRIANGLE,
new_or);
break;
case 4:
new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
new_v = faces[new_i]->GetVertices();
new_or = GetQuadOrientation(old_v, new_v);
dof_ord = fec->DofOrderForOrientation(Geometry::SQUARE,
new_or);
break;
}
old_dofs.SetSize(num_face_dofs);
new_dofs.SetSize(num_face_dofs);
for (j = 0; j < num_face_dofs; j++)
{
old_dofs[j] = offset + i * num_face_dofs + j;
new_dofs[j] = offset + new_i * num_face_dofs + dof_ord[j];
// we assumed the dofs are non-directional
// i.e. dof_ord[j] is >= 0
}
fes->DofsToVDofs(old_dofs);
fes->DofsToVDofs(new_dofs);
for (j = 0; j < old_dofs.Size(); j++)
(*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
}
offset += NumOfFaces * num_face_dofs;
delete faces_tbl;
}
// element dofs:
// element orientation may be different
if (num_elem_dofs > 0)
{
// matters when the 'fec' is
// - Pk on triangles, k >= 4
// - Qk on quads, k >= 3
// - Pk on tets, k >= 5
// - Qk on hexes, k >= 3
// - ...
}
// Update Tables, faces, etc
if (Dim > 2)
{
if (num_face_dofs == 0)
{
// needed for FE spaces that have face dofs, even if
// the 'Nodes' do not have face dofs.
GetElementToFaceTable();
GenerateFaces();
}
CheckBdrElementOrientation();
}
if (el_to_edge)
{
// update 'el_to_edge', 'be_to_edge' (2D), 'bel_to_edge' (3D)
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
if (Dim == 2)
{
// update 'faces' and 'faces_info'
GenerateFaces();
CheckBdrElementOrientation();
}
}
}
}
}
Mesh::Mesh(istream &input, Vector **data, int nprocessors,
int currentprocessor, int generate_edges)
{
Init();
InitTables();
Load(input, data, nprocessors, currentprocessor, generate_edges);
}
Mesh::Mesh(istream ** in, int np, int * dim)
{
int i,p,j,attr,ints[8];
int * shift = new int[np];
Dim = 3;
Init();
InitTables();
char first_line[1000];
for (p = 0; p < np; p++)
{
*in[p] >> first_line;
}
shift[0] = 0;
NumOfVertices = 0;
for (p = 0; p < np; p++)
{
*in[p] >> dim[p];
NumOfVertices += dim[p];
if (p < (np-1))
shift[p+1] = shift[p] + dim[p];
}
vertices.SetSize(NumOfVertices);
// Read the vertices
for (p = 0; p < np; p++)
{
for (i = 0; i < dim[p]; i++)
for (j = 0; j < Dim; j++)
*in[p] >> vertices[shift[p]+i](j);
}
int * edim = new int[np];
NumOfElements = 0;
for (p = 0; p < np; p++)
{
*in[p] >> edim[p];
NumOfElements += edim[p];
}
elements.SetSize(NumOfElements);
int e = 0;
// Read the elements
for (p = 0; p < np; p++)
for (i = 0; i < edim[p]; i++)
{
*in[p] >> attr;
for (j = 0; j<4; j++)
{
*in[p] >> ints[j];
ints[j]--;
ints[j] += shift[p];
}
#ifdef MFEM_USE_MEMALLOC
Tetrahedron *tet;
tet = TetMemory.Alloc();
tet->SetVertices(ints);
tet->SetAttribute(attr);
elements[e++] = tet;
#else
elements[e++] = new Tetrahedron(ints, attr);
#endif
}
int * bdim = new int[np];
NumOfBdrElements = 0;
for (p = 0; p < np; p++)
{
*in[p] >> bdim[p];
NumOfBdrElements += bdim[p];
}
boundary.SetSize(NumOfBdrElements);
int be=0;
// Read the boundary information.
for (p = 0; p < np; p++)
for (i = 0; i < bdim[p]; i++)
{
*in[p] >> attr;
for (j = 0; j < 3; j++)
{
*in[p] >> ints[j];
ints[j]--;
ints[j] += shift[p];
}
boundary[be++] = new Triangle(ints, attr);
}
SetAttributes();
delete [] shift;
delete [] bdim;
delete [] edim;
}
void Mesh::Load(istream &input, Vector **data, int nprocessors,
int currentprocessor, int generate_edges)
{
int i, j, ints[32], attr;
static int elems = 0, bdrelems = 0, totalNumberOfVertices = 0, np = 0;
#ifdef MFEM_DEBUG
if (!input)
mfem_error("Input file stream not opened : Mesh::Load");
#endif
Dim = 3;
InitTables();
int gid;
int totalNumberOfElements, totalNumberOfBdrElements;
int components, c;
if (data[1] == NULL)
components = 1;
else
components = 3;
// Read the vertices
input >> NumOfVertices >> totalNumberOfVertices;
if (np == 0)
{
vertices.SetSize(totalNumberOfVertices);
for (c = 0; c < components; c++)
data[c]->SetSize(totalNumberOfVertices);
}
for (i = 0; i < NumOfVertices; i++)
{
input >> gid;
for (j = 0; j < Dim; j++)
input >> vertices[gid](j);
for (c = 0; c < components; c++)
input >> (*(data[c]))(gid);
}
// Read the elements
input >> NumOfElements >> totalNumberOfElements;
if (np == 0)
elements.SetSize(totalNumberOfElements);
for (i = 0; i < NumOfElements; i++)
{
input >> attr;
for (j = 0; j < 4; j++)
input >> ints[j];
#ifdef MFEM_USE_MEMALLOC
Tetrahedron *tet;
tet = TetMemory.Alloc();
tet->SetVertices(ints);
tet->SetAttribute(attr);
elements[i+elems] = tet;
#else
elements[i+elems] = new Tetrahedron(ints, attr);
#endif
}
elems += NumOfElements;
// Read the boundary information.
input >> NumOfBdrElements >> totalNumberOfBdrElements;
if (np == 0)
boundary.SetSize(totalNumberOfBdrElements);
for (i = 0; i < NumOfBdrElements; i++)
{
input >> attr;
for (j = 0; j < 3; j++)
input >> ints[j];
boundary[i+bdrelems] = new Triangle(ints, attr);
}
bdrelems += NumOfBdrElements;
np++;
if (np == nprocessors)
{
NumOfVertices = totalNumberOfVertices;
NumOfElements = elems;
NumOfBdrElements = bdrelems;
MarkTetMeshForRefinement();
if (Dim == 3)
{
GetElementToFaceTable();
GenerateFaces();
CheckBdrElementOrientation();
}
if (generate_edges == 1)
{
el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
if (Dim == 2)
{
GenerateFaces();
CheckBdrElementOrientation();
}
}
else
NumOfEdges = 0;
SetAttributes();
np = elems = bdrelems = totalNumberOfVertices = 0;
}
}
void XYZ_VectorFunction(const Vector &p, Vector &v)
{
v = p;
}
void Mesh::SetNodalFESpace(FiniteElementSpace *nfes)
{
GridFunction *nodes = new GridFunction(nfes);
VectorFunctionCoefficient xyz(Dim, XYZ_VectorFunction);
nodes->ProjectCoefficient(xyz);
if (own_nodes) delete Nodes;
Nodes = nodes;
own_nodes = 1;
}
void Mesh::SetNodalGridFunction(GridFunction *nodes)
{
if (Nodes == NULL || Nodes->FESpace() != nodes->FESpace())
{
VectorFunctionCoefficient xyz(Dim, XYZ_VectorFunction);
nodes->ProjectCoefficient(xyz);
}
else
*nodes = *Nodes;
if (own_nodes) delete Nodes;
Nodes = nodes;
own_nodes = 0;
}
const FiniteElementSpace *Mesh::GetNodalFESpace()
{
return ((Nodes) ? Nodes->FESpace() : NULL);
}
void Mesh::CheckElementOrientation()
{
int i, j, k, wo = 0, *vi;
double *v[4];
if (Dim == 2)
{
DenseMatrix tri(2, 2);
for (i = 0; i < NumOfElements; i++)
{
vi = elements[i]->GetVertices();
for (j = 0; j < 3; j++)
v[j] = vertices[vi[j]]();
for (j = 0; j < 2; j++)
for (k = 0; k < 2; k++)
tri(j, k) = v[j+1][k] - v[0][k];
if (tri.Det() < 0.0)
switch (GetElementType(i))
{
case Element::TRIANGLE:
k = vi[0], vi[0] = vi[1], vi[1] = k, wo++;
break;
case Element::QUADRILATERAL:
k = vi[1], vi[1] = vi[3], vi[3] = k, wo++;
break;
}
}
}
if (Dim == 3)
{
DenseMatrix tet(3, 3);
for (i = 0; i < NumOfElements; i++)
{
vi = elements[i]->GetVertices();
switch (GetElementType(i))
{
case Element::TETRAHEDRON:
for (j = 0; j < 4; j++)
v[j] = vertices[vi[j]]();
for (j = 0; j < 3; j++)
for (k = 0; k < 3; k++)
tet(j, k) = v[j+1][k] - v[0][k];
if (tet.Det() < 0.0)
k = vi[0], vi[0] = vi[1], vi[1] = k, wo++;
break;
case Element::HEXAHEDRON:
// to do ...
break;
}
}
}
//#ifdef MFEM_DEBUG
if (wo > 0)
cout << "Orientation fixed in " << wo << " of "<< NumOfElements
<< " elements" << endl;
//#endif
}
int Mesh::GetTriOrientation(const int *base, const int *test)
{
int orient;
if (test[0] == base[0])
if (test[1] == base[1])
orient = 0; // (0, 1, 2)
else
orient = 5; // (0, 2, 1)
else if (test[0] == base[1])
if (test[1] == base[0])
orient = 1; // (1, 0, 2)
else
orient = 2; // (1, 2, 0)
else // test[0] == base[2]
if (test[1] == base[0])
orient = 4; // (2, 0, 1)
else
orient = 3; // (2, 1, 0)
#ifdef MFEM_DEBUG
static const int tri_orient[6][3] = {{0, 1, 2}, {1, 0, 2},
{2, 0, 1}, {2, 1, 0},
{1, 2, 0}, {0, 2, 1}};
const int *aor = tri_orient[orient];
for (int j = 0; j < 3; j++)
if (test[aor[j]] != base[j])
mfem_error("Mesh::GetTriOrientation(...)");
#endif
return orient;
}
int Mesh::GetQuadOrientation(const int *base, const int *test)
{
int i;
for (i = 0; i < 4; i++)
if (test[i] == base[0])
break;
#ifdef MFEM_DEBUG
static const int quad_orient[8][4] = {{0, 1, 2, 3}, {0, 3, 2, 1},
{1, 2, 3, 0}, {1, 0, 3, 2},
{2, 3, 0, 1}, {2, 1, 0, 3},
{3, 0, 1, 2}, {3, 2, 1, 0}};
int orient;
if (test[(i+1)%4] == base[1])
orient = 2*i;
else
orient = 2*i+1;
const int *aor = quad_orient[orient];
for (int j = 0; j < 4; j++)
if (test[aor[j]] != base[j])
{
cerr << "Mesh::GetQuadOrientation(...)" << endl;
cerr << " base = [";
for (int k = 0; k < 4; k++)
cerr << " " << base[k];
cerr << " ]\n test = [";
for (int k = 0; k < 4; k++)
cerr << " " << test[k];
cerr << " ]" << endl;
mfem_error();
}
#endif
if (test[(i+1)%4] == base[1])
return 2*i;
return 2*i+1;
}
void Mesh::CheckBdrElementOrientation()
{
int i, j, wo = 0;
if (Dim == 2)
{
for (i = 0; i < NumOfBdrElements; i++)
{
if (faces_info[be_to_edge[i]].Elem2No < 0) // boundary face
{
int *bv = boundary[i]->GetVertices();
int *fv = faces[be_to_edge[i]]->GetVertices();
if (bv[0] != fv[0])
{
j = bv[0]; bv[0] = bv[1]; bv[1] = j;
wo++;
}
}
}
}
if (Dim == 3)
{
int el, *bv, *ev;
int v[4];
for (i = 0; i < NumOfBdrElements; i++)
{
if (faces_info[be_to_face[i]].Elem2No < 0)
{ // boundary face
bv = boundary[i]->GetVertices();
el = faces_info[be_to_face[i]].Elem1No;
ev = elements[el]->GetVertices();
switch (GetElementType(el))
{
case Element::TETRAHEDRON:
{
int *fv = faces[be_to_face[i]]->GetVertices();
int orientation; // orientation of the bdr. elem. w.r.t. the
// corresponding face element (that's the base)
orientation = GetTriOrientation(fv, bv);
if (orientation % 2)
{
// wrong orientation -- swap vertices 0 and 1 so that
// we don't change the marked edge: (0,1,2) -> (1,0,2)
j = bv[0]; bv[0] = bv[1]; bv[1] = j;
wo++;
}
}
break;
case Element::HEXAHEDRON:
switch (faces_info[be_to_face[i]].Elem1Inf/64)
{
case 0:
v[0] = ev[3]; v[1] = ev[2]; v[2] = ev[1]; v[3] = ev[0];
break;
case 1:
v[0] = ev[0]; v[1] = ev[1]; v[2] = ev[5]; v[3] = ev[4];
break;
case 2:
v[0] = ev[1]; v[1] = ev[2]; v[2] = ev[6]; v[3] = ev[5];
break;
case 3:
v[0] = ev[2]; v[1] = ev[3]; v[2] = ev[7]; v[3] = ev[6];
break;
case 4:
v[0] = ev[3]; v[1] = ev[0]; v[2] = ev[4]; v[3] = ev[7];
break;
case 5:
v[0] = ev[4]; v[1] = ev[5]; v[2] = ev[6]; v[3] = ev[7];
break;
}
if (GetQuadOrientation(v, bv) % 2)
{
j = bv[0]; bv[0] = bv[2]; bv[2] = j;
wo++;
}
break;
}
}
}
}
//#ifdef MFEM_DEBUG
if (wo > 0)
cout << "Orientation fixed in " << wo << " of "<< NumOfBdrElements
<< " boundary elements" << endl;
//#endif
}
void Mesh::GetElementEdges(int i, Array<int> &edges, Array<int> &cor)
const
{
if (el_to_edge)
el_to_edge->GetRow(i, edges);
else
mfem_error("Mesh::GetElementEdges(...) element to edge table "
"is not generated.");
const int *v = elements[i]->GetVertices();
const int ne = elements[i]->GetNEdges();
cor.SetSize(ne);
for (int j = 0; j < ne; j++)
{
const int *e = elements[i]->GetEdgeVertices(j);
cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
}
}
void Mesh::GetBdrElementEdges(int i, Array<int> &edges, Array<int> &cor)
const
{
if (Dim == 2)
{
edges.SetSize(1);
cor.SetSize(1);
edges[0] = be_to_edge[i];
const int *v = boundary[i]->GetVertices();
cor[0] = (v[0]<v[1]) ? (1) : (-1);
}
else if (Dim == 3)
{
if (bel_to_edge)
bel_to_edge->GetRow(i, edges);
else
mfem_error("Mesh::GetBdrElementEdges(...)");
const int *v = boundary[i]->GetVertices();
const int ne = boundary[i]->GetNEdges();
cor.SetSize(ne);
for (int j = 0; j < ne; j++)
{
const int *e = boundary[i]->GetEdgeVertices(j);
cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
}
}
}
void Mesh::GetFaceEdges(int i, Array<int> &edges, Array<int> &o) const
{
if (Dim != 3)
return;
GetFaceEdgeTable(); // generate face_edge Table (if not generated)
face_edge->GetRow(i, edges);
// to do: set the orientation 'o'
}
void Mesh::GetEdgeVertices(int i, Array<int> &vert) const
{
if (Dim == 2 && faces.Size() == NumOfEdges)
{
faces[i]->GetVertices(vert);
}
else
{
GetEdgeVertexTable(); // generate edge_vertex Table (if not generated)
edge_vertex->GetRow(i, vert);
}
}
Table *Mesh::GetFaceEdgeTable() const
{
if (face_edge)
return face_edge;
if (Dim != 3)
return NULL;
DSTable v_to_v(NumOfVertices);
GetVertexToVertexTable(v_to_v);
face_edge = new Table;
GetElementArrayEdgeTable(faces, v_to_v, *face_edge);
return(face_edge);
}
Table *Mesh::GetEdgeVertexTable() const
{
if (edge_vertex)
return edge_vertex;
DSTable v_to_v(NumOfVertices);
GetVertexToVertexTable(v_to_v);
int nedges = v_to_v.NumberOfEntries();
edge_vertex = new Table(nedges, 2);
for (int i = 0; i < NumOfVertices; i++)
{
for (DSTable::RowIterator it(v_to_v, i); !it; ++it)
{
int j = it.Index();
edge_vertex->Push(j, i);
edge_vertex->Push(j, it.Column());
}
}
edge_vertex->Finalize();
return edge_vertex;
}
Table *Mesh::GetVertexToElementTable()
{
int i, j, nv, *v;
Table *vert_elem = new Table;
vert_elem->MakeI(NumOfVertices);
for (i = 0; i < NumOfElements; i++)
{
nv = elements[i]->GetNVertices();
v = elements[i]->GetVertices();
for (j = 0; j < nv; j++)
vert_elem->AddAColumnInRow(v[j]);
}
vert_elem->MakeJ();
for (i = 0; i < NumOfElements; i++)
{
nv = elements[i]->GetNVertices();
v = elements[i]->GetVertices();
for (j = 0; j < nv; j++)
vert_elem->AddConnection(v[j], i);
}
vert_elem->ShiftUpI();
return vert_elem;
}
void Mesh::GetElementFaces(int i, Array<int> &fcs, Array<int> &cor)
const
{
int n, j;
if (el_to_face)
el_to_face->GetRow(i, fcs);
else
mfem_error("Mesh::GetElementFaces(...) : el_to_face not generated.");
n = fcs.Size();
cor.SetSize(n);
for (j = 0; j < n; j++)
if (faces_info[fcs[j]].Elem1No == i)
cor[j] = faces_info[fcs[j]].Elem1Inf % 64;
#ifdef MFEM_DEBUG
else if (faces_info[fcs[j]].Elem2No == i)
cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
else
mfem_error("Mesh::GetElementFaces(...) : 2");
#else
else
cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
#endif
}
void Mesh::GetBdrElementFace(int i, int *f, int *o) const
{
const int *bv, *fv;
if (State == Mesh::TWO_LEVEL_COARSE)
{
// the coarse level 'be_to_face' and 'faces' are destroyed
mfem_error("Mesh::GetBdrElementFace (...)");
}
*f = be_to_face[i];
bv = boundary[i]->GetVertices();
fv = faces[be_to_face[i]]->GetVertices();
// find the orientation of the bdr. elem. w.r.t.
// the corresponding face element (that's the base)
switch (GetBdrElementType(i))
{
case Element::TRIANGLE:
*o = GetTriOrientation(fv, bv);
break;
case Element::QUADRILATERAL:
*o = GetQuadOrientation(fv, bv);
break;
default:
mfem_error("Mesh::GetBdrElementFace(...) 2");
}
}
int Mesh::GetFaceBaseGeometry(int i) const
{
// Here, we assume all faces are of the same type
switch (GetElementType(0))
{
case Element::TRIANGLE:
case Element::QUADRILATERAL:
return Geometry::SEGMENT; // in 2D 'face' is an edge
case Element::TETRAHEDRON:
return Geometry::TRIANGLE;
case Element::HEXAHEDRON:
return Geometry::SQUARE;
default:
mfem_error("Mesh::GetFaceBaseGeometry(...) #1");
}
return(-1);
#if 0
if (faces[i] == NULL)
switch (GetElementType(faces_info[i].Elem1No))
{
case Element::TETRAHEDRON:
return Geometry::TRIANGLE;
case Element::HEXAHEDRON:
return Geometry::SQUARE;
default:
mfem_error("Mesh::GetFaceBaseGeometry(...) #2");
}
else
return faces[i]->GetGeometryType();
#endif
}
int Mesh::GetBdrElementEdgeIndex(int i) const
{
if (Dim == 2)
return be_to_edge[i];
return be_to_face[i];
}
int Mesh::GetElementType(int i) const
{
Element *El = elements[i];
int t = El->GetType();
while (1)
if (t == Element::BISECTED ||
t == Element::QUADRISECTED ||
t == Element::OCTASECTED)
t = (El = ((RefinedElement *) El)->IAm())->GetType();
else
break;
return t;
}
int Mesh::GetBdrElementType(int i) const
{
Element *El = boundary[i];
int t = El->GetType();
while (1)
if (t == Element::BISECTED || t == Element::QUADRISECTED)
t = (El = ((RefinedElement *) El)->IAm())->GetType();
else
break;
return t;
}
void Mesh::GetPointMatrix(int i, DenseMatrix &pointmat) const
{
int k, j, nv;
const int *v;
v = elements[i]->GetVertices();
nv = elements[i]->GetNVertices();
pointmat.SetSize(Dim, nv);
for (k = 0; k < Dim; k++)
for (j = 0; j < nv; j++)
pointmat(k, j) = vertices[v[j]](k);
}
void Mesh::GetBdrPointMatrix(int i,DenseMatrix &pointmat) const
{
int k, j, nv;
const int *v;
v = boundary[i]->GetVertices();
nv = boundary[i]->GetNVertices();
pointmat.SetSize(Dim, nv);
for (k = 0; k < Dim; k++)
for (j = 0; j < nv; j++)
pointmat(k, j) = vertices[v[j]](k);
}
double Mesh::GetLength(int i, int j) const
{
const double *vi = vertices[i]();
const double *vj = vertices[j]();
double length = 0.;
for (int k = 0; k < Dim; k++)
length += (vi[k]-vj[k])*(vi[k]-vj[k]);
return sqrt(length);
}
// static method
void Mesh::GetElementArrayEdgeTable(const Array<Element*> &elem_array,
const DSTable &v_to_v, Table &el_to_edge)
{
el_to_edge.MakeI(elem_array.Size());
for (int i = 0; i < elem_array.Size(); i++)
{
el_to_edge.AddColumnsInRow(i, elem_array[i]->GetNEdges());
}
el_to_edge.MakeJ();
for (int i = 0; i < elem_array.Size(); i++)
{
const int *v = elem_array[i]->GetVertices();
const int ne = elem_array[i]->GetNEdges();
for (int j = 0; j < ne; j++)
{
const int *e = elem_array[i]->GetEdgeVertices(j);
el_to_edge.AddConnection(i, v_to_v(v[e[0]], v[e[1]]));
}
}
el_to_edge.ShiftUpI();
}
void Mesh::GetVertexToVertexTable(DSTable & v_to_v) const
{
for (int i = 0; i < NumOfElements; i++)
{
const int *v = elements[i]->GetVertices();
const int ne = elements[i]->GetNEdges();
for (int j = 0; j < ne; j++)
{
const int *e = elements[i]->GetEdgeVertices(j);
v_to_v.Push(v[e[0]], v[e[1]]);
}
}
}
int Mesh::GetElementToEdgeTable(Table & e_to_f, Array<int> &be_to_f)
{
int i, NumberOfEdges;
DSTable v_to_v(NumOfVertices);
GetVertexToVertexTable(v_to_v);
NumberOfEdges = v_to_v.NumberOfEntries();
// Fill the element to edge table
GetElementArrayEdgeTable(elements, v_to_v, e_to_f);
if (Dim == 2)
{
// Initialize the indeces for the boundary elements.
be_to_f.SetSize(NumOfBdrElements);
for (i = 0; i < NumOfBdrElements; i++)
{
const int *v = boundary[i]->GetVertices();
be_to_f[i] = v_to_v(v[0], v[1]);
}
}
else if (Dim == 3)
{
if (bel_to_edge == NULL)
bel_to_edge = new Table;
GetElementArrayEdgeTable(boundary, v_to_v, *bel_to_edge);
}
else
mfem_error("1D GetElementToEdgeTable is not yet implemented.");
// Return the number of edges
return NumberOfEdges;
}
const Table & Mesh::ElementToElementTable()
{
if (el_to_el)
return *el_to_el;
if (Dim == 2)
{
Table edge_el;
Transpose(ElementToEdgeTable(), edge_el);
el_to_el = new Table(NumOfElements, 4); // 4 is the max. # of edges
for (int i = 0; i < edge_el.Size(); i++)
if (edge_el.RowSize(i) > 1)
{
const int *el = edge_el.GetRow(i);
el_to_el->Push(el[0], el[1]);
el_to_el->Push(el[1], el[0]);
}
el_to_el->Finalize();
}
else if (Dim == 3)
{
el_to_el = new Table(NumOfElements, 6); // 6 is the max. # of faces
for (int i = 0; i < faces_info.Size(); i++)
if (faces_info[i].Elem2No >= 0)
{
el_to_el->Push(faces_info[i].Elem1No, faces_info[i].Elem2No);
el_to_el->Push(faces_info[i].Elem2No, faces_info[i].Elem1No);
}
el_to_el->Finalize();
}
else
mfem_error("Mesh::ElementToElementTable() in 1D is not implemented!");
return *el_to_el;
}
const Table & Mesh::ElementToFaceTable() const
{
if (el_to_face == NULL)
mfem_error("Mesh::ElementToFaceTable()");
return *el_to_face;
}
const Table & Mesh::ElementToEdgeTable() const
{
if (el_to_edge == NULL)
mfem_error("Mesh::ElementToEdgeTable()");
return *el_to_edge;
}
void Mesh::AddSegmentFaceElement(int lf, int gf, int el, int v0, int v1)
{
if (faces[gf] == NULL) // this will be elem1
{
faces[gf] = new Segment(v0, v1);
faces_info[gf].Elem1No = el;
faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
faces_info[gf].Elem2No = -1; // in case there's no other side
}
else // this will be elem2
{
#ifdef MFEM_DEBUG
int *v = faces[gf]->GetVertices();
if (v[1] != v0 || v[0] != v1)
mfem_error("Mesh::AddSegmentFaceElement(...)");
#endif
faces_info[gf].Elem2No = el;
faces_info[gf].Elem2Inf = 64 * lf + 1;
}
}
void Mesh::AddTriangleFaceElement(int lf, int gf, int el,
int v0, int v1, int v2)
{
if (faces[gf] == NULL) // this will be elem1
{
faces[gf] = new Triangle(v0, v1, v2);
faces_info[gf].Elem1No = el;
faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
faces_info[gf].Elem2No = -1; // in case there's no other side
}
else // this will be elem2
{
int orientation, vv[3] = { v0, v1, v2 };
orientation = GetTriOrientation(faces[gf]->GetVertices(), vv);
#ifdef MFEM_DEBUG
if (orientation % 2 == 0)
mfem_error("Mesh::AddTriangleFaceElement(...)");
#endif
faces_info[gf].Elem2No = el;
faces_info[gf].Elem2Inf = 64 * lf + orientation;
}
}
void Mesh::AddQuadFaceElement(int lf, int gf, int el,
int v0, int v1, int v2, int v3)
{
if (faces_info[gf].Elem1No < 0) // this will be elem1
{
faces[gf] = new Quadrilateral(v0, v1, v2, v3);
faces_info[gf].Elem1No = el;
faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
faces_info[gf].Elem2No = -1; // in case there's no other side
}
else // this will be elem2
{
int vv[4] = { v0, v1, v2, v3 };
int oo = GetQuadOrientation(faces[gf]->GetVertices(), vv);
#ifdef MFEM_DEBUG
if (oo % 2 == 0)
mfem_error("Mesh::AddQuadFaceElement(...)");
#endif
faces_info[gf].Elem2No = el;
faces_info[gf].Elem2Inf = 64 * lf + oo;
}
}
void Mesh::GenerateFaces()
{
int i, nfaces;
nfaces = (Dim == 2) ? NumOfEdges : NumOfFaces;
for (i = 0; i < faces.Size(); i++)
FreeElement(faces[i]);
// (re)generate the interior faces and the info for them
faces.SetSize(nfaces);
faces_info.SetSize(nfaces);
for (i = 0; i < nfaces; i++)
{
faces[i] = NULL;
faces_info[i].Elem1No = -1;
}
for (i = 0; i < NumOfElements; i++)
{
const int *v = elements[i]->GetVertices();
const int *ef;
if (Dim == 2)
{
ef = el_to_edge->GetRow(i);
const int ne = elements[i]->GetNEdges();
for (int j = 0; j < ne; j++)
{
const int *e = elements[i]->GetEdgeVertices(j);
AddSegmentFaceElement(j, ef[j], i, v[e[1]], v[e[0]]);
}
}
else
{
ef = el_to_face->GetRow(i);
switch (GetElementType(i))
{
case Element::TETRAHEDRON:
AddTriangleFaceElement(0, ef[0], i, v[1], v[2], v[3]);
AddTriangleFaceElement(1, ef[1], i, v[0], v[3], v[2]);
AddTriangleFaceElement(2, ef[2], i, v[0], v[1], v[3]);
AddTriangleFaceElement(3, ef[3], i, v[0], v[2], v[1]);
break;
case Element::HEXAHEDRON:
AddQuadFaceElement(0, ef[0], i, v[3], v[2], v[1], v[0]);
AddQuadFaceElement(1, ef[1], i, v[0], v[1], v[5], v[4]);
AddQuadFaceElement(2, ef[2], i, v[1], v[2], v[6], v[5]);
AddQuadFaceElement(3, ef[3], i, v[2], v[3], v[7], v[6]);
AddQuadFaceElement(4, ef[4], i, v[3], v[0], v[4], v[7]);
AddQuadFaceElement(5, ef[5], i, v[4], v[5], v[6], v[7]);
break;
}
}
}
}
STable3D *Mesh::GetElementToFaceTable(int ret_ftbl)
{
int i, *v;
STable3D *faces_tbl;
if (el_to_face != NULL)
delete el_to_face;
el_to_face = new Table(NumOfElements, 6); // must be 6 for hexahedra
faces_tbl = new STable3D(NumOfVertices);
for (i = 0; i < NumOfElements; i++)
{
v = elements[i]->GetVertices();
switch (GetElementType(i))
{
case Element::TETRAHEDRON:
el_to_face->Push(i, faces_tbl->Push(v[1], v[2], v[3]));
el_to_face->Push(i, faces_tbl->Push(v[0], v[3], v[2]));
el_to_face->Push(i, faces_tbl->Push(v[0], v[1], v[3]));
el_to_face->Push(i, faces_tbl->Push(v[0], v[2], v[1]));
break;
case Element::HEXAHEDRON:
// find the face by the vertices with the smallest 3 numbers
// z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
el_to_face->Push(i, faces_tbl->Push4(v[3], v[2], v[1], v[0]));
el_to_face->Push(i, faces_tbl->Push4(v[0], v[1], v[5], v[4]));
el_to_face->Push(i, faces_tbl->Push4(v[1], v[2], v[6], v[5]));
el_to_face->Push(i, faces_tbl->Push4(v[2], v[3], v[7], v[6]));
el_to_face->Push(i, faces_tbl->Push4(v[3], v[0], v[4], v[7]));
el_to_face->Push(i, faces_tbl->Push4(v[4], v[5], v[6], v[7]));
break;
}
}
el_to_face->Finalize();
NumOfFaces = faces_tbl->NumberOfElements();
be_to_face.SetSize(NumOfBdrElements);
for (i = 0; i < NumOfBdrElements; i++)
{
v = boundary[i]->GetVertices();
switch (GetBdrElementType(i))
{
case Element::TRIANGLE:
be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
break;
case Element::QUADRILATERAL:
be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
break;
}
}
if (ret_ftbl)
return faces_tbl;
delete faces_tbl;
return NULL;
}
#ifdef MFEM_USE_MPI
// auxiliary function for qsort
static int mfem_less(const void *x, const void *y)
{
if (*(int*)x < *(int*)y)
return 1;
if (*(int*)x > *(int*)y)
return -1;
return 0;
}
// METIS prototypes
typedef int idxtype;
extern "C" {
void METIS_PartGraphRecursive(int*, idxtype*, idxtype*, idxtype*, idxtype*,
int*, int*, int*, int*, int*, idxtype*);
void METIS_PartGraphKway(int*, idxtype*, idxtype*, idxtype*, idxtype*,
int*, int*, int*, int*, int*, idxtype*);
void METIS_PartGraphVKway(int*, idxtype*, idxtype*, idxtype*, idxtype*,
int*, int*, int*, int*, int*, idxtype*);
}
#endif
int *Mesh::GeneratePartitioning(int nparts, int part_method)
{
#ifdef MFEM_USE_MPI
int i, *partitioning;
ElementToElementTable();
partitioning = new int[NumOfElements];
if (nparts == 1)
{
for (i = 0; i < NumOfElements; i++)
partitioning[i] = 0;
}
else
{
int *I, *J, n;
int wgtflag = 0;
int numflag = 0;
int options[5];
int edgecut;
n = NumOfElements;
I = el_to_el->GetI();
J = el_to_el->GetJ();
options[0] = 0;
// Sort the neighbor lists
if (part_method >= 0 && part_method <= 2)
for (i = 0; i < n; i++)
qsort(&J[I[i]], I[i+1]-I[i], sizeof(int), &mfem_less);
// This function should be used to partition a graph into a small
// number of partitions (less than 8).
if (part_method == 0 || part_method == 3)
METIS_PartGraphRecursive(&n,
(idxtype *) I,
(idxtype *) J,
(idxtype *) NULL,
(idxtype *) NULL,
&wgtflag,
&numflag,
&nparts,
options,
&edgecut,
(idxtype *) partitioning);
// This function should be used to partition a graph into a large
// number of partitions (greater than 8).
if (part_method == 1 || part_method == 4)
METIS_PartGraphKway(&n,
(idxtype *) I,
(idxtype *) J,
(idxtype *) NULL,
(idxtype *) NULL,
&wgtflag,
&numflag,
&nparts,
options,
&edgecut,
(idxtype *) partitioning);
// The objective of this partitioning is to minimize the total
// communication volume
if (part_method == 2 || part_method == 5)
METIS_PartGraphVKway(&n,
(idxtype *) I,
(idxtype *) J,
(idxtype *) NULL,
(idxtype *) NULL,
&wgtflag,
&numflag,
&nparts,
options,
&edgecut,
(idxtype *) partitioning);
#ifdef MFEM_DEBUG
cout << "Mesh::GeneratePartitioning(...): edgecut = "
<< edgecut << endl;
#endif
}
if (el_to_el)
delete el_to_el;
el_to_el = NULL;
// Check for empty partitionings (a "feature" in METIS)
{
Array< Pair<int,int> > psize(nparts);
for (i = 0; i < nparts; i++)
{
psize[i].one = 0;
psize[i].two = i;
}
for (i = 0; i < NumOfElements; i++)
psize[partitioning[i]].one++;
int empty_parts = 0;
for (i = 0; i < nparts; i++)
if (psize[i].one == 0)
empty_parts++;
// This code just split the largest partitionings in two.
// Do we need to replace it with something better?
if (empty_parts)
{
cerr << "Mesh::GeneratePartitioning returned " << empty_parts
<< " empty parts!" << endl;
SortPairs<int,int>(psize, nparts);
for (i = nparts-1; i > nparts-1-empty_parts; i--)
psize[i].one /= 2;
for (int j = 0; j < NumOfElements; j++)
for (i = nparts-1; i > nparts-1-empty_parts; i--)
if (psize[i].one == 0 || partitioning[j] != psize[i].two)
continue;
else
{
partitioning[j] = psize[nparts-1-i].two;
psize[i].one--;
}
}
}
return partitioning;
#else
mfem_error("Mesh::GeneratePartitioning(...): "
"MFEM was compiled without Metis.");
return NULL;
#endif
}
/* required: 0 <= partitioning[i] < num_part */
void FindPartitioningComponents(Table &elem_elem,
const Array<int> &partitioning,
Array<int> &component,
Array<int> &num_comp)
{
int i, j, k;
int num_elem, *i_elem_elem, *j_elem_elem;
num_elem = elem_elem.Size();
i_elem_elem = elem_elem.GetI();
j_elem_elem = elem_elem.GetJ();
component.SetSize(num_elem);
Array<int> elem_stack(num_elem);
int stack_p, stack_top_p, elem;
int num_part;
num_part = -1;
for (i = 0; i < num_elem; i++)
{
if (partitioning[i] > num_part)
num_part = partitioning[i];
component[i] = -1;
}
num_part++;
num_comp.SetSize(num_part);
for (i = 0; i < num_part; i++)
num_comp[i] = 0;
stack_p = 0;
stack_top_p = 0; // points to the first unused element in the stack
for (elem = 0; elem < num_elem; elem++)
{
if (component[elem] >= 0)
continue;
component[elem] = num_comp[partitioning[elem]]++;
elem_stack[stack_top_p++] = elem;
for ( ; stack_p < stack_top_p; stack_p++)
{
i = elem_stack[stack_p];
for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
{
k = j_elem_elem[j];
if (partitioning[k] == partitioning[i])
if (component[k] < 0)
{
component[k] = component[i];
elem_stack[stack_top_p++] = k;
}
else if (component[k] != component[i])
{
mfem_error("FindPartitioningComponents");
}
}
}
}
}
void Mesh::CheckPartitioning(int *partitioning)
{
int i, n_empty, n_mcomp;
Array<int> component, num_comp;
const Array<int> _partitioning(partitioning, GetNE());
ElementToElementTable();
FindPartitioningComponents(*el_to_el, _partitioning, component, num_comp);
n_empty = n_mcomp = 0;
for (i = 0; i < num_comp.Size(); i++)
if (num_comp[i] == 0)
n_empty++;
else if (num_comp[i] > 1)
n_mcomp++;
if (n_empty > 0)
{
cout << "Mesh::CheckPartitioning(...) :\n"
<< "The following subdomains are empty :\n";
for (i = 0; i < num_comp.Size(); i++)
if (num_comp[i] == 0)
cout << ' ' << i;
cout << endl;
}
if (n_mcomp > 0)
{
cout << "Mesh::CheckPartitioning(...) :\n"
<< "The following subdomains are NOT connected :\n";
for (i = 0; i < num_comp.Size(); i++)
if (num_comp[i] > 1)
cout << ' ' << i;
cout << endl;
}
if (n_empty == 0 && n_mcomp == 0)
cout << "Mesh::CheckPartitioning(...) : "
"All subdomains are connected." << endl;
if (el_to_el)
delete el_to_el;
el_to_el = NULL;
}
// compute the coefficients of the polynomial in t:
// c(0)+c(1)*t+...+c(d)*t^d = det(A+t*B)
// where A, B are (d x d), d=2,3
void DetOfLinComb(const DenseMatrix &A, const DenseMatrix &B, Vector &c)
{
const double *a = A.Data();
const double *b = B.Data();
c.SetSize(A.Size()+1);
switch (A.Size())
{
case 2:
{
// det(A+t*B) = |a0 a2| / |a0 b2| + |b0 a2| \
// |a1 a3| + \ |a1 b3| |b1 a3| / * t +
// |b0 b2|
// |b1 b3| * t^2
c(0) = a[0]*a[3]-a[1]*a[2];
c(1) = a[0]*b[3]-a[1]*b[2]+b[0]*a[3]-b[1]*a[2];
c(2) = b[0]*b[3]-b[1]*b[2];
}
break;
case 3:
{
// |a0 a3 a6|
// det(A+t*B) = |a1 a4 a7| +
// |a2 a5 a8|
// / |b0 a3 a6| |a0 b3 a6| |a0 a3 b6| \
// | |b1 a4 a7| + |a1 b4 a7| + |a1 a4 b7| | * t +
// \ |b2 a5 a8| |a2 b5 a8| |a2 a5 b8| /
// / |a0 b3 b6| |b0 a3 b6| |b0 b3 a6| \
// | |a1 b4 b7| + |b1 a4 b7| + |b1 b4 a7| | * t^2 +
// \ |a2 b5 b8| |b2 a5 b8| |b2 b5 a8| /
// |b0 b3 b6|
// |b1 b4 b7| * t^3
// |b2 b5 b8|
c(0) = (a[0] * (a[4] * a[8] - a[5] * a[7]) +
a[1] * (a[5] * a[6] - a[3] * a[8]) +
a[2] * (a[3] * a[7] - a[4] * a[6]));
c(1) = (b[0] * (a[4] * a[8] - a[5] * a[7]) +
b[1] * (a[5] * a[6] - a[3] * a[8]) +
b[2] * (a[3] * a[7] - a[4] * a[6]) +
a[0] * (b[4] * a[8] - b[5] * a[7]) +
a[1] * (b[5] * a[6] - b[3] * a[8]) +
a[2] * (b[3] * a[7] - b[4] * a[6]) +
a[0] * (a[4] * b[8] - a[5] * b[7]) +
a[1] * (a[5] * b[6] - a[3] * b[8]) +
a[2] * (a[3] * b[7] - a[4] * b[6]));
c(2) = (a[0] * (b[4] * b[8] - b[5] * b[7]) +
a[1] * (b[5] * b[6] - b[3] * b[8]) +
a[2] * (b[3] * b[7] - b[4] * b[6]) +
b[0] * (a[4] * b[8] - a[5] * b[7]) +
b[1] * (a[5] * b[6] - a[3] * b[8]) +
b[2] * (a[3] * b[7] - a[4] * b[6]) +
b[0] * (b[4] * a[8] - b[5] * a[7]) +
b[1] * (b[5] * a[6] - b[3] * a[8]) +
b[2] * (b[3] * a[7] - b[4] * a[6]));
c(3) = (b[0] * (b[4] * b[8] - b[5] * b[7]) +
b[1] * (b[5] * b[6] - b[3] * b[8]) +
b[2] * (b[3] * b[7] - b[4] * b[6]));
}
break;
default:
mfem_error("DetOfLinComb(...)");
}
}
inline void swap(double &a, double &b)
{
const double t = a;
a = b;
b = t;
}
// compute the real roots of
// z(0)+z(1)*x+...+z(d)*x^d = 0, d=2,3;
// the roots are returned in x, sorted in increasing order;
// it is assumed that x is at least if size d;
// return the number of roots counting multiplicity;
// return -1 if all z(i) are 0.
int FindRoots(const Vector &z, Vector &x)
{
int d = z.Size()-1;
if (d > 3 || d < 0)
mfem_error("FindRoots(...)");
while (z(d) == 0.0)
{
if (d == 0)
return(-1);
d--;
}
switch (d)
{
case 0:
{
return 0;
}
case 1:
{
x(0) = -z(0)/z(1);
return 1;
}
case 2:
{
double a = z(2), b = z(1), c = z(0);
double D = b*b-4*a*c;
if (D < 0.0)
{
return 0;
}
if (D == 0.0)
{
x(0) = x(1) = -0.5 * b / a;
return 2; // root with multiplicity 2
}
if (b == 0.0)
{
x(0) = -(x(1) = fabs(0.5 * sqrt(D) / a));
return 2;
}
else
{
double t;
if (b > 0.0)
t = -0.5 * (b + sqrt(D));
else
t = -0.5 * (b - sqrt(D));
x(0) = t / a;
x(1) = c / t;
if (x(0) > x(1))
swap(x(0), x(1));
return 2;
}
}
case 3:
{
double a = z(2)/z(3), b = z(1)/z(3), c = z(0)/z(3);
// find the real roots of x^3 + a x^2 + b x + c = 0
double Q = (a * a - 3 * b) / 9;
double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
double Q3 = Q * Q * Q;
double R2 = R * R;
if (R2 == Q3)
{
if (Q == 0)
{
x(0) = x(1) = x(2) = - a / 3;
}
else
{
double sqrtQ = sqrt(Q);
if (R > 0)
{
x(0) = -2 * sqrtQ - a / 3;
x(1) = x(2) = sqrtQ - a / 3;
}
else
{
x(0) = x(1) = - sqrtQ - a / 3;
x(2) = 2 * sqrtQ - a / 3;
}
}
return 3;
}
else if (R2 < Q3)
{
double theta = acos(R / sqrt(Q3));
double A = -2 * sqrt(Q);
double x0, x1, x2;
x0 = A * cos(theta / 3) - a / 3;
x1 = A * cos((theta + 2.0 * M_PI) / 3) - a / 3;
x2 = A * cos((theta - 2.0 * M_PI) / 3) - a / 3;
/* Sort x0, x1, x2 */
if (x0 > x1)
swap(x0, x1);
if (x1 > x2)
{
swap(x1, x2);
if (x0 > x1)
swap(x0, x1);
}
x(0) = x0;
x(1) = x1;
x(2) = x2;
return 3;
}
else
{
double A;
if (R >= 0.0)
A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
else
A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
x(0) = A + Q / A - a / 3;
return 1;
}
}
}
return 0;
}
void FindTMax(Vector &c, Vector &x, double &tmax,
const double factor, const int Dim)
{
const double c0 = c(0);
c(0) = c0 * (1.0 - pow(factor, -Dim));
int nr = FindRoots(c, x);
for (int j = 0; j < nr; j++)
{
if (x(j) > tmax)
break;
if (x(j) >= 0.0)
{
tmax = x(j);
break;
}
}
c(0) = c0 * (1.0 - pow(factor, Dim));
nr = FindRoots(c, x);
for (int j = 0; j < nr; j++)
{
if (x(j) > tmax)
break;
if (x(j) >= 0.0)
{
tmax = x(j);
break;
}
}
}
void Mesh::CheckDisplacements(const Vector &displacements, double &tmax)
{
int nvs = vertices.Size();
DenseMatrix P, V, DS, PDS(Dim), VDS(Dim);
Vector c(Dim+1), x(Dim);
const double factor = 2.0;
// check for tangling assuming constant speed
if (tmax < 1.0)
tmax = 1.0;
for (int i = 0; i < NumOfElements; i++)
{
Element *el = elements[i];
int nv = el->GetNVertices();
int *v = el->GetVertices();
P.SetSize(Dim, nv);
V.SetSize(Dim, nv);
for (int j = 0; j < Dim; j++)
for (int k = 0; k < nv; k++)
{
P(j, k) = vertices[v[k]](j);
V(j, k) = displacements(v[k]+j*nvs);
}
DS.SetSize(nv, Dim);
const FiniteElement *fe =
GetTransformationFEforElementType(el->GetType());
// check if det(P.DShape+t*V.DShape) > 0 for all x and 0<=t<=1
switch (el->GetType())
{
case Element::TRIANGLE:
case Element::TETRAHEDRON:
{
// DS is constant
fe->CalcDShape(Geometries.GetCenter(fe->GetGeomType()), DS);
Mult(P, DS, PDS);
Mult(V, DS, VDS);
DetOfLinComb(PDS, VDS, c);
if (c(0) <= 0.0)
tmax = 0.0;
else
FindTMax(c, x, tmax, factor, Dim);
}
break;
case Element::QUADRILATERAL:
{
const IntegrationRule &ir = fe->GetNodes();
for (int j = 0; j < nv; j++)
{
fe->CalcDShape(ir.IntPoint(j), DS);
Mult(P, DS, PDS);
Mult(V, DS, VDS);
DetOfLinComb(PDS, VDS, c);
if (c(0) <= 0.0)
tmax = 0.0;
else
FindTMax(c, x, tmax, factor, Dim);
}
}
break;
default:
mfem_error("Mesh::CheckDisplacements(...)");
}
}
}
void Mesh::MoveVertices(const Vector &displacements)
{
for (int i = 0, nv = vertices.Size(); i < nv; i++)
for (int j = 0; j < Dim; j++)
vertices[i](j) += displacements(j*nv+i);
}
void Mesh::GetVertices(Vector &vert_coord) const
{
int nv = vertices.Size();
vert_coord.SetSize(nv*Dim);
for (int i = 0; i < nv; i++)
for (int j = 0; j < Dim; j++)
vert_coord(j*nv+i) = vertices[i](j);
}
void Mesh::SetVertices(const Vector &vert_coord)
{
for (int i = 0, nv = vertices.Size(); i < nv; i++)
for (int j = 0; j < Dim; j++)
vertices[i](j) = vert_coord(j*nv+i);
}
void Mesh::MoveNodes(const Vector &displacements)
{
if (Nodes)
(*Nodes) += displacements;
else
MoveVertices(displacements);
}
void Mesh::GetNodes(Vector &node_coord) const
{
if (Nodes)
node_coord = (*Nodes);
else
GetVertices(node_coord);
}
void Mesh::SetNodes(const Vector &node_coord)
{
if (Nodes)
(*Nodes) = node_coord;
else
SetVertices(node_coord);
}
void Mesh::NewNodes(GridFunction &nodes)
{
if (own_nodes) delete Nodes;
Nodes = &nodes;
own_nodes = 0;
}
void Mesh::AverageVertices(int * indexes, int n, int result)
{
int j, k;
for (k = 0; k < Dim; k++)
vertices[result](k) = vertices[indexes[0]](k);
for (j = 1; j < n; j++)
for (k = 0; k < Dim; k++)
vertices[result](k) += vertices[indexes[j]](k);
for (k = 0; k < Dim; k++)
vertices[result](k) *= (1.0 / n);
}
void Mesh::UpdateNodes()
{
FiniteElementSpace *cfes = Nodes->FESpace()->SaveUpdate();
SparseMatrix *R =
Nodes->FESpace()->GlobalRestrictionMatrix(cfes, 0);
delete cfes;
{
Vector cNodes = *Nodes;
Nodes->Update();
R->MultTranspose(cNodes, *Nodes);
}
delete R;
SetState(Mesh::TWO_LEVEL_FINE);
// update the vertices?
}
void Mesh::QuadUniformRefinement()
{
int i, j, *v, vv[2], attr, wtls = WantTwoLevelState;
const int *e;
if (Nodes) // curved mesh
{
UseTwoLevelState(1);
}
SetState(Mesh::NORMAL);
if (el_to_edge == NULL)
{
el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
}
int oedge = NumOfVertices;
int oelem = oedge + NumOfEdges;
DeleteCoarseTables();
vertices.SetSize(oelem + NumOfElements);
for (i = 0; i < NumOfElements; i++)
{
v = elements[i]->GetVertices();
AverageVertices(v, 4, oelem+i);
e = el_to_edge->GetRow(i);
vv[0] = v[0], vv[1] = v[1]; AverageVertices(vv, 2, oedge+e[0]);
vv[0] = v[1], vv[1] = v[2]; AverageVertices(vv, 2, oedge+e[1]);
vv[0] = v[2], vv[1] = v[3]; AverageVertices(vv, 2, oedge+e[2]);
vv[0] = v[3], vv[1] = v[0]; AverageVertices(vv, 2, oedge+e[3]);
}
elements.SetSize(4 * NumOfElements);
for (i = 0; i < NumOfElements; i++)
{
attr = elements[i]->GetAttribute();
v = elements[i]->GetVertices();
e = el_to_edge->GetRow(i);
j = NumOfElements + 3 * i;
elements[j+0] = new Quadrilateral(oedge+e[0], v[1], oedge+e[1],
oelem+i, attr);
elements[j+1] = new Quadrilateral(oelem+i, oedge+e[1], v[2],
oedge+e[2], attr);
elements[j+2] = new Quadrilateral(oedge+e[3], oelem+i, oedge+e[2],
v[3], attr);
if (WantTwoLevelState)
{
QuadrisectedElement *qe;
qe = new QuadrisectedElement(elements[i]->Duplicate());
qe->FirstChild = elements[i];
qe->Child2 = j;
qe->Child3 = j+1;
qe->Child4 = j+2;
elements[i] = qe;
}
v[1] = oedge+e[0];
v[2] = oelem+i;
v[3] = oedge+e[3];
}
boundary.SetSize(2 * NumOfBdrElements);
for (i = 0; i < NumOfBdrElements; i++)
{
attr = boundary[i]->GetAttribute();
v = boundary[i]->GetVertices();
j = NumOfBdrElements + i;
boundary[j] = new Segment(oedge+be_to_edge[i], v[1], attr);
if (WantTwoLevelState)
{
#ifdef MFEM_USE_MEMALLOC
BisectedElement *be = BEMemory.Alloc();
#else
BisectedElement *be = new BisectedElement;
#endif
be->SetCoarseElem(boundary[i]->Duplicate());
be->FirstChild = boundary[i];
be->SecondChild = j;
boundary[i] = be;
}
v[1] = oedge+be_to_edge[i];
}
if (WantTwoLevelState)
{
c_NumOfVertices = NumOfVertices;
c_NumOfEdges = NumOfEdges;
c_NumOfElements = NumOfElements;
c_NumOfBdrElements = NumOfBdrElements;
RefinedElement::State = RefinedElement::FINE;
State = Mesh::TWO_LEVEL_FINE;
}
NumOfVertices = oelem + NumOfElements;
NumOfElements = 4 * NumOfElements;
NumOfBdrElements = 2 * NumOfBdrElements;
NumOfFaces = 0;
if (WantTwoLevelState)
{
f_NumOfVertices = NumOfVertices;
f_NumOfElements = NumOfElements;
f_NumOfBdrElements = NumOfBdrElements;
}
if (el_to_edge != NULL)
{
if (WantTwoLevelState)
{
c_el_to_edge = el_to_edge;
Swap(be_to_edge, fc_be_to_edge); // save coarse be_to_edge
f_el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
el_to_edge = f_el_to_edge;
f_NumOfEdges = NumOfEdges;
}
else
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
GenerateFaces();
}
#ifdef MFEM_DEBUG
CheckElementOrientation();
CheckBdrElementOrientation();
#endif
if (Nodes) // curved mesh
{
UpdateNodes();
UseTwoLevelState(wtls);
}
}
void Mesh::HexUniformRefinement()
{
int i, wtls = WantTwoLevelState;
int * v;
const int *e, *f;
int vv[4];
if (Nodes) // curved mesh
{
UseTwoLevelState(1);
}
SetState(Mesh::NORMAL);
if (el_to_edge == NULL)
{
el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
}
if (el_to_face == NULL)
GetElementToFaceTable();
int oedge = NumOfVertices;
int oface = oedge + NumOfEdges;
int oelem = oface + NumOfFaces;
DeleteCoarseTables();
vertices.SetSize(oelem + NumOfElements);
for (i = 0; i < NumOfElements; i++)
{
v = elements[i]->GetVertices();
AverageVertices(v, 8, oelem+i);
f = el_to_face->GetRow(i);
vv[0] = v[3], vv[1] = v[2], vv[2] = v[1], vv[3] = v[0];
AverageVertices(vv, 4, oface+f[0]);
vv[0] = v[0], vv[1] = v[1], vv[2] = v[5], vv[3] = v[4];
AverageVertices(vv, 4, oface+f[1]);
vv[0] = v[1], vv[1] = v[2], vv[2] = v[6], vv[3] = v[5];
AverageVertices(vv, 4, oface+f[2]);
vv[0] = v[2], vv[1] = v[3], vv[2] = v[7], vv[3] = v[6];
AverageVertices(vv, 4, oface+f[3]);
vv[0] = v[3], vv[1] = v[0], vv[2] = v[4], vv[3] = v[7];
AverageVertices(vv, 4, oface+f[4]);
vv[0] = v[4], vv[1] = v[5], vv[2] = v[6], vv[3] = v[7];
AverageVertices(vv, 4, oface+f[5]);
e = el_to_edge->GetRow(i);
vv[0] = v[0], vv[1] = v[1]; AverageVertices(vv, 2, oedge+e[0]);
vv[0] = v[1], vv[1] = v[2]; AverageVertices(vv, 2, oedge+e[1]);
vv[0] = v[2], vv[1] = v[3]; AverageVertices(vv, 2, oedge+e[2]);
vv[0] = v[3], vv[1] = v[0]; AverageVertices(vv, 2, oedge+e[3]);
vv[0] = v[4], vv[1] = v[5]; AverageVertices(vv, 2, oedge+e[4]);
vv[0] = v[5], vv[1] = v[6]; AverageVertices(vv, 2, oedge+e[5]);
vv[0] = v[6], vv[1] = v[7]; AverageVertices(vv, 2, oedge+e[6]);
vv[0] = v[7], vv[1] = v[4]; AverageVertices(vv, 2, oedge+e[7]);
vv[0] = v[0], vv[1] = v[4]; AverageVertices(vv, 2, oedge+e[8]);
vv[0] = v[1], vv[1] = v[5]; AverageVertices(vv, 2, oedge+e[9]);
vv[0] = v[2], vv[1] = v[6]; AverageVertices(vv, 2, oedge+e[10]);
vv[0] = v[3], vv[1] = v[7]; AverageVertices(vv, 2, oedge+e[11]);
}
int attr, j, k;
elements.SetSize(8 * NumOfElements);
for (i = 0; i < NumOfElements; i++)
{
attr = elements[i]->GetAttribute();
v = elements[i]->GetVertices();
e = el_to_edge->GetRow(i);
f = el_to_face->GetRow(i);
j = NumOfElements + 7 * i;
elements[j+0] = new Hexahedron(oedge+e[0], v[1], oedge+e[1], oface+f[0],
oface+f[1], oedge+e[9], oface+f[2],
oelem+i, attr);
elements[j+1] = new Hexahedron(oface+f[0], oedge+e[1], v[2], oedge+e[2],
oelem+i, oface+f[2], oedge+e[10],
oface+f[3], attr);
elements[j+2] = new Hexahedron(oedge+e[3], oface+f[0], oedge+e[2], v[3],
oface+f[4], oelem+i, oface+f[3],
oedge+e[11], attr);
elements[j+3] = new Hexahedron(oedge+e[8], oface+f[1], oelem+i,
oface+f[4], v[4], oedge+e[4], oface+f[5],
oedge+e[7], attr);
elements[j+4] = new Hexahedron(oface+f[1], oedge+e[9], oface+f[2],
oelem+i, oedge+e[4], v[5], oedge+e[5],
oface+f[5], attr);
elements[j+5] = new Hexahedron(oelem+i, oface+f[2], oedge+e[10],
oface+f[3], oface+f[5], oedge+e[5], v[6],
oedge+e[6], attr);
elements[j+6] = new Hexahedron(oface+f[4], oelem+i, oface+f[3],
oedge+e[11], oedge+e[7], oface+f[5],
oedge+e[6], v[7], attr);
if (WantTwoLevelState)
{
OctasectedElement *oe;
oe = new OctasectedElement(elements[i]->Duplicate());
oe->FirstChild = elements[i];
for (k = 0; k < 7; k++)
oe->Child[k] = j + k;
elements[i] = oe;
}
v[1] = oedge+e[0];
v[2] = oface+f[0];
v[3] = oedge+e[3];
v[4] = oedge+e[8];
v[5] = oface+f[1];
v[6] = oelem+i;
v[7] = oface+f[4];
}
boundary.SetSize(4 * NumOfBdrElements);
for (i = 0; i < NumOfBdrElements; i++)
{
attr = boundary[i]->GetAttribute();
v = boundary[i]->GetVertices();
e = bel_to_edge->GetRow(i);
f = & be_to_face[i];
j = NumOfBdrElements + 3 * i;
boundary[j+0] = new Quadrilateral(oedge+e[0], v[1], oedge+e[1],
oface+f[0], attr);
boundary[j+1] = new Quadrilateral(oface+f[0], oedge+e[1], v[2],
oedge+e[2], attr);
boundary[j+2] = new Quadrilateral(oedge+e[3], oface+f[0], oedge+e[2],
v[3], attr);
if (WantTwoLevelState)
{
QuadrisectedElement *qe;
qe = new QuadrisectedElement(boundary[i]->Duplicate());
qe->FirstChild = boundary[i];
qe->Child2 = j;
qe->Child3 = j+1;
qe->Child4 = j+2;
boundary[i] = qe;
}
v[1] = oedge+e[0];
v[2] = oface+f[0];
v[3] = oedge+e[3];
}
if (WantTwoLevelState)
{
c_NumOfVertices = NumOfVertices;
c_NumOfEdges = NumOfEdges;
c_NumOfFaces = NumOfFaces;
c_NumOfElements = NumOfElements;
c_NumOfBdrElements = NumOfBdrElements;
RefinedElement::State = RefinedElement::FINE;
State = Mesh::TWO_LEVEL_FINE;
}
NumOfVertices = oelem + NumOfElements;
NumOfElements = 8 * NumOfElements;
NumOfBdrElements = 4 * NumOfBdrElements;
if (WantTwoLevelState)
{
f_NumOfVertices = NumOfVertices;
f_NumOfElements = NumOfElements;
f_NumOfBdrElements = NumOfBdrElements;
}
if (el_to_edge != NULL)
{
if (WantTwoLevelState)
{
c_el_to_edge = el_to_edge;
f_el_to_edge = new Table;
c_bel_to_edge = bel_to_edge;
bel_to_edge = NULL;
NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
el_to_edge = f_el_to_edge;
f_bel_to_edge = bel_to_edge;
f_NumOfEdges = NumOfEdges;
}
else
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
}
if (el_to_face != NULL)
{
if (WantTwoLevelState)
{
c_el_to_face = el_to_face;
el_to_face = NULL;
Swap(faces_info, fc_faces_info);
}
GetElementToFaceTable();
GenerateFaces();
if (WantTwoLevelState)
{
f_el_to_face = el_to_face;
f_NumOfFaces = NumOfFaces;
}
}
#ifdef MFEM_DEBUG
CheckBdrElementOrientation();
#endif
if (Nodes) // curved mesh
{
UpdateNodes();
UseTwoLevelState(wtls);
}
// When 'WantTwoLevelState' is true the coarse level
// 'be_to_face' and 'faces'
// are destroyed !!!
}
void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
{
int i, j, ind, nedges, wtls = WantTwoLevelState;
Array<int> v;
if (Nodes) // curved mesh
{
UseTwoLevelState(1);
}
SetState(Mesh::NORMAL);
DeleteCoarseTables();
if (Dim == 1) // --------------------------------------------------------
{
int cne = NumOfElements, cnv = NumOfVertices;
NumOfVertices += marked_el.Size();
NumOfElements += marked_el.Size();
vertices.SetSize(NumOfVertices);
elements.SetSize(NumOfElements);
for (j = 0; j < marked_el.Size(); j++)
{
i = marked_el[j];
int *vert = elements[i]->GetVertices();
vertices[cnv+j](0) = 0.5 * ( vertices[vert[0]](0) +
vertices[vert[1]](0) );
elements[cne+j] = new Segment(cnv+j, vert[1],
elements[i]->GetAttribute());
vert[1] = cnv+j;
}
} // end of 'if (Dim == 1)'
else if (Dim == 2) // ---------------------------------------------------
{
if (WantTwoLevelState)
{
c_NumOfVertices = NumOfVertices;
c_NumOfEdges = NumOfEdges;
c_NumOfElements = NumOfElements;
c_NumOfBdrElements = NumOfBdrElements;
}
// 1. Get table of vertex to vertex connections.
DSTable v_to_v(NumOfVertices);
GetVertexToVertexTable(v_to_v);
// 2. Get edge to element connections in arrays edge1 and edge2
nedges = v_to_v.NumberOfEntries();
int *edge1 = new int[nedges];
int *edge2 = new int[nedges];
int *middle = new int[nedges];
for (i = 0; i < nedges; i++)
edge1[i] = edge2[i] = middle[i] = -1;
for (i = 0; i < NumOfElements; i++)
{
elements[i]->GetVertices(v);
for (j = 1; j < v.Size(); j++)
{
ind = v_to_v(v[j-1], v[j]);
(edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
}
ind = v_to_v(v[0], v[v.Size()-1]);
(edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
}
// 3. Do the red refinement.
for (i = 0; i < marked_el.Size(); i++)
RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
// 4. Do the green refinement (to get conforming mesh).
int need_refinement;
do
{
need_refinement = 0;
for (i = 0; i < nedges; i++)
if (middle[i] != -1 && edge1[i] != -1)
{
need_refinement = 1;
GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
}
}
while (need_refinement == 1);
// 5. Update the boundary elements.
int v1[2], v2[2], bisect, temp;
temp = NumOfBdrElements;
for (i = 0; i < temp; i++)
{
boundary[i]->GetVertices(v);
bisect = v_to_v(v[0], v[1]);
if (middle[bisect] != -1) // the element was refined (needs updating)
{
if (boundary[i]->GetType() == Element::SEGMENT)
{
v1[0] = v[0]; v1[1] = middle[bisect];
v2[0] = middle[bisect]; v2[1] = v[1];
if (WantTwoLevelState)
{
boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
#ifdef MFEM_USE_MEMALLOC
BisectedElement *aux = BEMemory.Alloc();
aux->SetCoarseElem(boundary[i]);
#else
BisectedElement *aux = new BisectedElement(boundary[i]);
#endif
aux->FirstChild =
new Segment(v1, boundary[i]->GetAttribute());
aux->SecondChild = NumOfBdrElements;
boundary[i] = aux;
NumOfBdrElements++;
}
else
{
boundary[i]->SetVertices(v1);
boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
}
}
else
mfem_error("Only bisection of segment is implemented"
" for bdr elem.");
}
}
NumOfBdrElements = boundary.Size();
// 6. Free the allocated memory.
delete [] edge1;
delete [] edge2;
delete [] middle;
#ifdef MFEM_DEBUG
CheckElementOrientation();
#endif
if (WantTwoLevelState)
{
f_NumOfVertices = NumOfVertices;
f_NumOfElements = NumOfElements;
f_NumOfBdrElements = NumOfBdrElements;
RefinedElement::State = RefinedElement::FINE;
State = Mesh::TWO_LEVEL_FINE;
}
if (el_to_edge != NULL)
{
if (WantTwoLevelState)
{
c_el_to_edge = el_to_edge;
Swap(be_to_edge, fc_be_to_edge); // save coarse be_to_edge
f_el_to_edge = new Table;
NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
el_to_edge = f_el_to_edge;
f_NumOfEdges = NumOfEdges;
}
else
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
GenerateFaces();
}
}
else if (Dim == 3) // ---------------------------------------------------
{
if (WantTwoLevelState)
{
c_NumOfVertices = NumOfVertices;
c_NumOfEdges = NumOfEdges;
c_NumOfFaces = NumOfFaces;
c_NumOfElements = NumOfElements;
c_NumOfBdrElements = NumOfBdrElements;
}
// 1. Get table of vertex to vertex connections.
DSTable v_to_v(NumOfVertices);
GetVertexToVertexTable(v_to_v);
// 2. Get edge to element connections in arrays edge1 and edge2
nedges = v_to_v.NumberOfEntries();
int *middle = new int[nedges];
for (i = 0; i < nedges; i++)
middle[i] = -1;
// 3. Do the red refinement.
int ii;
switch (type)
{
case 1:
for (i = 0; i < marked_el.Size(); i++)
Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
break;
case 2:
for (i = 0; i < marked_el.Size(); i++)
{
Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
}
break;
case 3:
for (i = 0; i < marked_el.Size(); i++)
{
Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
ii = NumOfElements - 1;
Bisection(ii, v_to_v, NULL, NULL, middle);
Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
Bisection(ii, v_to_v, NULL, NULL, middle);
Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
Bisection(NumOfElements-1, v_to_v, NULL, NULL, middle);
Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
}
break;
}
if (WantTwoLevelState)
{
RefinedElement::State = RefinedElement::FINE;
State = Mesh::TWO_LEVEL_FINE;
}
// 4. Do the green refinement (to get conforming mesh).
int need_refinement;
// int need_refinement, onoe, max_gen = 0;
do
{
// int redges[2], type, flag;
need_refinement = 0;
// onoe = NumOfElements;
// for (i = 0; i < onoe; i++)
for (i = 0; i < NumOfElements; i++)
{
// ((Tetrahedron *)elements[i])->ParseRefinementFlag(redges, type, flag);
// if (flag > max_gen) max_gen = flag;
if (elements[i]->NeedRefinement(v_to_v, middle))
{
need_refinement = 1;
Bisection(i, v_to_v, NULL, NULL, middle);
}
}
}
while (need_refinement == 1);
// cout << "Maximum generation: " << max_gen << endl;
// 5. Update the boundary elements.
do
{
need_refinement = 0;
for (i = 0; i < NumOfBdrElements; i++)
if (boundary[i]->NeedRefinement(v_to_v, middle))
{
need_refinement = 1;
Bisection(i, v_to_v, middle);
}
}
while (need_refinement == 1);
// 6. Un-mark the Pf elements.
int refinement_edges[2], type, flag;
for (i = 0; i < NumOfElements; i++)
{
Element *El = elements[i];
while (El->GetType() == Element::BISECTED)
El = ((BisectedElement *) El)->FirstChild;
((Tetrahedron *)El)->ParseRefinementFlag(refinement_edges, type,
flag);
if (type == Tetrahedron::TYPE_PF)
((Tetrahedron *)El)->CreateRefinementFlag(refinement_edges,
Tetrahedron::TYPE_PU,
flag);
}
NumOfBdrElements = boundary.Size();
// 7. Free the allocated memory.
delete [] middle;
#ifdef MFEM_DEBUG
CheckElementOrientation();
#endif
if (el_to_edge != NULL)
{
if (WantTwoLevelState)
{
c_el_to_edge = el_to_edge;
f_el_to_edge = new Table;
c_bel_to_edge = bel_to_edge;
bel_to_edge = NULL;
NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
el_to_edge = f_el_to_edge;
f_bel_to_edge = bel_to_edge;
}
else
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
}
if (el_to_face != NULL)
{
if (WantTwoLevelState)
{
c_el_to_face = el_to_face;
el_to_face = NULL;
Swap(faces_info, fc_faces_info);
}
GetElementToFaceTable();
GenerateFaces();
if (WantTwoLevelState)
{
f_el_to_face = el_to_face;
}
}
if (WantTwoLevelState)
{
f_NumOfVertices = NumOfVertices;
f_NumOfEdges = NumOfEdges;
f_NumOfFaces = NumOfFaces;
f_NumOfElements = NumOfElements;
f_NumOfBdrElements = NumOfBdrElements;
}
} // end 'if (Dim == 3)'
if (Nodes) // curved mesh
{
UpdateNodes();
UseTwoLevelState(wtls);
}
}
void Mesh::UniformRefinement()
{
if (meshgen == 1)
{
Array<int> elem_to_refine(GetNE());
for (int i = 0; i < elem_to_refine.Size(); i++)
elem_to_refine[i] = i;
LocalRefinement(elem_to_refine);
}
else if (Dim == 2)
QuadUniformRefinement();
else if (Dim == 3)
HexUniformRefinement();
else
mfem_error("Mesh::UniformRefinement()");
}
void Mesh::Bisection(int i, const DSTable &v_to_v,
int *edge1, int *edge2, int *middle)
{
int *vert;
int v[2][4], v_new, bisect, t;
Element **pce;
Vertex V;
if (WantTwoLevelState)
{
pce = &(elements[i]);
while (1)
{
t = pce[0]->GetType();
if (t == Element::BISECTED)
pce = & ( ((BisectedElement *) pce[0])->FirstChild );
else if (t == Element::QUADRISECTED)
pce = & ( ((QuadrisectedElement *) pce[0])->FirstChild );
else
break;
}
}
else
t = elements[i]->GetType();
if (t == Element::TRIANGLE)
{
Triangle *tri;
if (WantTwoLevelState)
tri = (Triangle *) pce[0];
else
tri = (Triangle *) elements[i];
vert = tri->GetVertices();
// 1. Get the index for the new vertex in v_new.
bisect = v_to_v(vert[0], vert[1]);
#ifdef MFEM_DEBUG
if (bisect < 0)
mfem_error("Mesh::Bisection(...): ERROR");
#endif
if (middle[bisect] == -1)
{
v_new = NumOfVertices++;
V(0) = 0.5 * (vertices[vert[0]](0) + vertices[vert[1]](0));
V(1) = 0.5 * (vertices[vert[0]](1) + vertices[vert[1]](1));
V(2) = 0.0;
vertices.Append(V);
// Put the element that may need refinement (because of this
// bisection) in edge1, or -1 if no more refinement is needed.
if (edge1[bisect] == i)
edge1[bisect] = edge2[bisect];
middle[bisect] = v_new;
}
else
{
v_new = middle[bisect];
// This edge will require no more refinement.
edge1[bisect] = -1;
}
// 2. Set the node indices for the new elements in v[0] and v[1] so that
// the edge marked for refinement is between the first two nodes.
v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
if (WantTwoLevelState)
{
#ifdef MFEM_USE_MEMALLOC
BisectedElement *aux = BEMemory.Alloc();
aux->SetCoarseElem(tri);
#else
BisectedElement *aux = new BisectedElement(tri);
#endif
aux->FirstChild = tri = new Triangle(v[0], tri->GetAttribute());
aux->SecondChild = NumOfElements;
pce[0] = aux;
}
else
tri->SetVertices(v[0]); // changes vert[0..2] !!!
elements.Append(new Triangle(v[1], tri->GetAttribute()));
// 3. edge1 and edge2 may have to be changed for the second triangle.
if (v[1][0] < v_to_v.NumberOfRows() && v[1][1] < v_to_v.NumberOfRows())
{
bisect = v_to_v(v[1][0], v[1][1]);
#ifdef MFEM_DEBUG
if (bisect < 0)
mfem_error("Mesh::Bisection(...): ERROR 2");
#endif
if (edge1[bisect] == i)
edge1[bisect] = NumOfElements;
else if (edge2[bisect] == i)
edge2[bisect] = NumOfElements;
}
NumOfElements++;
}
else if (t == Element::TETRAHEDRON)
{
int j, type, new_type, old_redges[2], new_redges[2][2], flag;
Tetrahedron *tet;
if (WantTwoLevelState)
tet = (Tetrahedron *) pce[0];
else
tet = (Tetrahedron *) elements[i];
if (tet->GetRefinementFlag() == 0)
mfem_error("Mesh::Bisection : TETRAHEDRON element is not marked for "
"refinement.");
vert = tet->GetVertices();
// 1. Get the index for the new vertex in v_new.
bisect = v_to_v(vert[0], vert[1]);
if (bisect == -1)
{
tet->ParseRefinementFlag(old_redges, type, flag);
cerr << "Error in Bisection(...) of tetrahedron!" << endl
<< " redge[0] = " << old_redges[0]
<< " redge[1] = " << old_redges[1]
<< " type = " << type
<< " flag = " << flag << endl;
mfem_error();
}
if (middle[bisect] == -1)
{
v_new = NumOfVertices++;
for (j = 0; j < 3; j++)
V(j) = 0.5 * (vertices[vert[0]](j) + vertices[vert[1]](j));
vertices.Append(V);
middle[bisect] = v_new;
}
else
v_new = middle[bisect];
// 2. Set the node indices for the new elements in v[2][4] so that
// the edge marked for refinement is between the first two nodes.
tet->ParseRefinementFlag(old_redges, type, flag);
v[0][3] = v_new;
v[1][3] = v_new;
new_redges[0][0] = 2;
new_redges[0][1] = 1;
new_redges[1][0] = 2;
new_redges[1][1] = 1;
switch (old_redges[0])
{
case 2:
v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
if (type == Tetrahedron::TYPE_PF) new_redges[0][1] = 4;
break;
case 3:
v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
break;
case 5:
v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
}
switch (old_redges[1])
{
case 1:
v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
if (type == Tetrahedron::TYPE_PF) new_redges[1][0] = 3;
break;
case 4:
v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
break;
case 5:
v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
}
int attr = tet->GetAttribute();
if (WantTwoLevelState)
{
#ifdef MFEM_USE_MEMALLOC
BisectedElement *aux = BEMemory.Alloc();
aux->SetCoarseElem(tet);
tet = TetMemory.Alloc();
tet->SetVertices(v[0]);
tet->SetAttribute(attr);
#else
BisectedElement *aux = new BisectedElement(tet);
tet = new Tetrahedron(v[0], attr);
#endif
aux->FirstChild = tet;
aux->SecondChild = NumOfElements;
pce[0] = aux;
}
else
tet->SetVertices(v[0]);
// 'tet' now points to the first child
{
#ifdef MFEM_USE_MEMALLOC
Tetrahedron *tet2 = TetMemory.Alloc();
tet2->SetVertices(v[1]);
tet2->SetAttribute(attr);
elements.Append(tet2);
#else
elements.Append(new Tetrahedron(v[1], attr));
#endif
}
// 3. Set the bisection flag
switch (type)
{
case Tetrahedron::TYPE_PU:
new_type = Tetrahedron::TYPE_PF; break;
case Tetrahedron::TYPE_PF:
new_type = Tetrahedron::TYPE_A; break;
default:
new_type = Tetrahedron::TYPE_PU;
}
tet->CreateRefinementFlag(new_redges[0], new_type, flag+1);
((Tetrahedron *)elements[NumOfElements])->
CreateRefinementFlag(new_redges[1], new_type, flag+1);
NumOfElements++;
}
else
mfem_error("Bisection for now works only for triangles & tetrahedra.");
}
void Mesh::Bisection(int i, const DSTable &v_to_v, int *middle)
{
int *vert;
int v[2][3], v_new, bisect, t;
Element **pce;
if (WantTwoLevelState)
{
pce = &(boundary[i]);
while (1)
{
t = pce[0]->GetType();
if (t == Element::BISECTED)
pce = & ( ((BisectedElement *) pce[0])->FirstChild );
else if (t == Element::QUADRISECTED)
pce = & ( ((QuadrisectedElement *) pce[0])->FirstChild );
else
break;
}
}
else
t = boundary[i]->GetType();
if (t == Element::TRIANGLE)
{
Triangle *tri;
if (WantTwoLevelState)
tri = (Triangle *) pce[0];
else
tri = (Triangle *) boundary[i];
vert = tri->GetVertices();
// 1. Get the index for the new vertex in v_new.
bisect = v_to_v(vert[0], vert[1]);
if (middle[bisect] == -1)
mfem_error("Error in Bisection(...) of boundary triangle!");
else
v_new = middle[bisect];
// 2. Set the node indices for the new elements in v[0] and v[1] so that
// the edge marked for refinement is between the first two nodes.
v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
if (WantTwoLevelState)
{
#ifdef MFEM_USE_MEMALLOC
BisectedElement *aux = BEMemory.Alloc();
aux->SetCoarseElem(tri);
#else
BisectedElement *aux = new BisectedElement(tri);
#endif
aux->FirstChild = tri = new Triangle(v[0], tri->GetAttribute());
aux->SecondChild = NumOfBdrElements;
pce[0] = aux;
}
else
boundary[i]->SetVertices(v[0]);
// 'tri' now points to the first child
boundary.Append(new Triangle(v[1], tri->GetAttribute()));
NumOfBdrElements++;
}
else
mfem_error("Bisection of boundary elements works only for triangles!");
}
void Mesh::UniformRefinement(int i, const DSTable &v_to_v,
int *edge1, int *edge2, int *middle)
{
Array<int> v;
int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
double coord[2];
if (elements[i]->GetType() == Element::TRIANGLE)
{
elements[i]->GetVertices(v);
// 1. Get the indeces for the new vertices in array v_new
bisect[0] = v_to_v(v[0],v[1]);
bisect[1] = v_to_v(v[1],v[2]);
bisect[2] = v_to_v(v[0],v[2]);
#ifdef MFEM_DEBUG
if (bisect[0] < 0 || bisect[1] < 0 || bisect[2] < 0)
mfem_error("Mesh::UniformRefinement(...): ERROR");
#endif
for (j = 0; j < 3; j++) // for the 3 edges fix v_new
if (middle[bisect[j]] == -1)
{
v_new[j] = NumOfVertices++;
coord[0] = (vertices[v[j]](0) + vertices[v[(j+1)%3]](0))/2.;
coord[1] = (vertices[v[j]](1) + vertices[v[(j+1)%3]](1))/2.;
Vertex V(coord[0], coord[1]);
vertices.Append(V);
// Put the element that may need refinement (because of this
// bisection) in edge1, or -1 if no more refinement is needed.
if (edge1[bisect[j]] == i)
edge1[bisect[j]] = edge2[bisect[j]];
middle[bisect[j]] = v_new[j];
}
else
{
v_new[j] = middle[bisect[j]];
// This edge will require no more refinement.
edge1[bisect[j]] = -1;
}
// 2. Set the node indeces for the new elements in v1, v2, v3 & v4 so that
// the edges marked for refinement be between the first two nodes.
v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
elements.Append(new Triangle(v1, elements[i]->GetAttribute()));
elements.Append(new Triangle(v2, elements[i]->GetAttribute()));
elements.Append(new Triangle(v3, elements[i]->GetAttribute()));
if (WantTwoLevelState)
{
QuadrisectedElement *aux = new QuadrisectedElement(elements[i]);
aux->FirstChild = new Triangle(v4, elements[i]->GetAttribute());
aux->Child2 = NumOfElements;
aux->Child3 = NumOfElements+1;
aux->Child4 = NumOfElements+2;
elements[i] = aux;
}
else
{
elements[i]->SetVertices(v4);
}
NumOfElements += 3;
}
else
mfem_error("Uniform refinement for now works only for triangles.");
}
void Mesh::SetState(int s)
{
if (State != Mesh::NORMAL && s == Mesh::NORMAL)
{
// two level state --->> normal state
int i, t;
for (i = 0; i < f_NumOfElements; )
{
t = elements[i]->GetType();
if (t == Element::BISECTED ||
t == Element::QUADRISECTED ||
t == Element::OCTASECTED)
{
RefinedElement *aux = (RefinedElement *) elements[i];
elements[i] = aux->FirstChild;
FreeElement(aux->CoarseElem);
FreeElement(aux);
}
else
i++;
}
for (i = 0; i < f_NumOfBdrElements; )
{
t = boundary[i]->GetType();
if (t == Element::BISECTED ||
t == Element::QUADRISECTED ||
t == Element::OCTASECTED)
{
RefinedElement *aux = (RefinedElement *) boundary[i];
boundary[i] = aux->FirstChild;
FreeElement(aux->CoarseElem);
FreeElement(aux);
}
else
i++;
}
if (el_to_edge != NULL)
{
delete c_el_to_edge;
el_to_edge = f_el_to_edge;
if (Dim == 2)
{
if (State == Mesh::TWO_LEVEL_COARSE)
Swap(be_to_edge, fc_be_to_edge);
fc_be_to_edge.DeleteAll();
}
if (Dim == 3)
{
delete c_bel_to_edge;
bel_to_edge = f_bel_to_edge;
}
}
if (el_to_face != NULL)
{
delete c_el_to_face;
el_to_face = f_el_to_face;
if (State == Mesh::TWO_LEVEL_COARSE)
Swap(faces_info, fc_faces_info);
fc_faces_info.DeleteAll();
}
NumOfVertices = f_NumOfVertices;
NumOfEdges = f_NumOfEdges;
NumOfFaces = f_NumOfFaces;
NumOfElements = f_NumOfElements;
NumOfBdrElements = f_NumOfBdrElements;
RefinedElement::State = RefinedElement::FINE;
State = s;
}
else if (State == Mesh::TWO_LEVEL_COARSE && s == Mesh::TWO_LEVEL_FINE)
{
if (el_to_edge != NULL)
{
el_to_edge = f_el_to_edge;
if (Dim == 2)
Swap(be_to_edge, fc_be_to_edge);
if (Dim == 3)
bel_to_edge = f_bel_to_edge;
}
if (el_to_face != NULL)
{
el_to_face = f_el_to_face;
Swap(faces_info, fc_faces_info);
}
NumOfVertices = f_NumOfVertices;
NumOfEdges = f_NumOfEdges;
NumOfFaces = f_NumOfFaces;
NumOfElements = f_NumOfElements;
NumOfBdrElements = f_NumOfBdrElements;
RefinedElement::State = RefinedElement::FINE;
State = s;
}
else if (State == Mesh::TWO_LEVEL_FINE && s == Mesh::TWO_LEVEL_COARSE)
{
if (el_to_edge != NULL)
{
el_to_edge = c_el_to_edge;
if (Dim == 2)
Swap(be_to_edge, fc_be_to_edge);
if (Dim == 3)
bel_to_edge = c_bel_to_edge;
}
if (el_to_face != NULL)
{
el_to_face = c_el_to_face;
Swap(faces_info, fc_faces_info);
}
NumOfVertices = c_NumOfVertices;
NumOfEdges = c_NumOfEdges;
NumOfFaces = c_NumOfFaces;
NumOfElements = c_NumOfElements;
NumOfBdrElements = c_NumOfBdrElements;
RefinedElement::State = RefinedElement::COARSE;
State = s;
}
else if (State != s)
mfem_error("Oops! Mesh::SetState");
}
int Mesh::GetNumFineElems(int i)
{
int t;
if (Dim == 2)
{
t = elements[i]->GetType();
if (t == Element::QUADRISECTED)
return 4;
else if (t == Element::BISECTED)
{
// assuming that the elements are either BisectedElements or
// regular elements
int n = 1;
BisectedElement *aux = (BisectedElement *) elements[i];
do
{
n += GetNumFineElems(aux->SecondChild);
if (aux->FirstChild->GetType() != Element::BISECTED)
break;
aux = (BisectedElement *) (aux->FirstChild);
}
while (1);
return n;
}
}
else if (Dim == 3)
{
// assuming that the element is a BisectedElement,
// OctasectedElement (with children that are regular elements) or
// regular element
t = elements[i]->GetType();
if (t == Element::BISECTED)
{
int n = 1;
BisectedElement *aux = (BisectedElement *) elements[i];
do
{
n += GetNumFineElems (aux->SecondChild);
if (aux->FirstChild->GetType() != Element::BISECTED)
break;
aux = (BisectedElement *) (aux->FirstChild);
}
while (1);
return n;
}
else if (t == Element::OCTASECTED)
return 8;
return 1; // regular element (i.e. it is not refined)
}
return 1; // the element is not refined
}
int Mesh::GetBisectionHierarchy(Element *E)
{
if (E->GetType() == Element::BISECTED)
{
int L, R, n, s, lb, rb;
L = GetBisectionHierarchy(((BisectedElement *)E)->FirstChild);
n = ((BisectedElement *)E)->SecondChild;
R = GetBisectionHierarchy(elements[n]);
n = 1; s = 1;
lb = rb = 1;
do
{
int nlb, nrb;
nlb = nrb = 0;
while (lb > 0)
{
n |= ((L & 1) << s);
s++;
nlb += (L & 1);
L = (L >> 1);
lb--;
}
while (rb > 0)
{
n |= ((R & 1) << s);
s++;
nrb += (R & 1);
R = (R >> 1);
rb--;
}
lb = 2 * nlb; rb = 2 * nrb;
}
while (lb > 0 || rb > 0);
return n;
}
return 0;
}
int Mesh::GetRefinementType(int i)
{
int t;
if (Dim == 2)
{
t = elements[i]->GetType();
if (t == Element::QUADRISECTED)
{
t = ((QuadrisectedElement *)elements[i])->CoarseElem->GetType();
if (t == Element::QUADRILATERAL)
return 1; // refinement type for quadrisected QUADRILATERAL
else
return 2; // refinement type for quadrisected TRIANGLE
}
else if (t == Element::BISECTED)
{
int type;
type = GetBisectionHierarchy(elements[i]);
if (type == 0)
mfem_error("Mesh::GetRefinementType(...)");
return type+2;
}
}
else if (Dim == 3)
{
int redges[2], type, flag;
Element *E = elements[i];
Tetrahedron *tet;
t = E->GetType();
if (t != Element::BISECTED)
if (t == Element::OCTASECTED)
return 1; // refinement type for octasected CUBE
else
return 0;
// Bisected TETRAHEDRON
tet = (Tetrahedron *) (((BisectedElement *) E)->CoarseElem);
tet->ParseRefinementFlag(redges, type, flag);
if (type == Tetrahedron::TYPE_A && redges[0] == 2)
type = 5;
else if (type == Tetrahedron::TYPE_M && redges[0] == 2)
type = 6;
type++;
type |= ( GetBisectionHierarchy(E) << 3 );
if (type < 8) type = 0;
return type;
}
return 0; // no refinement
}
int Mesh::GetFineElem(int i, int j)
{
int t;
if (Dim == 2)
{
t = elements[i]->GetType();
if (t == Element::QUADRISECTED)
{
QuadrisectedElement *aux = (QuadrisectedElement *) elements[i];
if (aux->CoarseElem->GetType() == Element::QUADRILATERAL)
switch (j)
{
case 0: return i;
case 1: return aux->Child2;
case 2: return aux->Child3;
case 3: return aux->Child4;
default: *((int *)NULL) = -1; // crash it
}
else // quadrisected TRIANGLE
switch (j)
{
case 0: return aux->Child2;
case 1: return aux->Child3;
case 2: return aux->Child4;
case 3: return i;
default: *((int *)NULL) = -1; // crash it
}
}
else if (t == Element::BISECTED)
{
int n = 0;
BisectedElement *aux = (BisectedElement *) elements[i];
do
{
int k = GetFineElem(aux->SecondChild, j-n);
if (k >= 0)
return k;
n -= k; // (-k) is the number of the leaves in this SecondChild
// n is the number of the leaves in
// the SecondChild-ren so far
if (aux->FirstChild->GetType() != Element::BISECTED)
break;
aux = (BisectedElement *) (aux->FirstChild);
}
while (1);
if (j > n) // i.e. if (j >= n+1)
return -(n+1);
return i; // j == n, i.e. j is the index of the last leaf
}
}
else if (Dim == 3)
{
t = elements[i]->GetType();
if (t == Element::BISECTED)
{
int n = 0;
BisectedElement *aux = (BisectedElement *) elements[i];
do
{
int k = GetFineElem(aux->SecondChild, j-n);
if (k >= 0)
return k;
n -= k; // (-k) is the number of the leaves in this SecondChild
// n is the number of the leaves in
// the SecondChild-ren so far
if (aux->FirstChild->GetType() != Element::BISECTED)
break;
aux = (BisectedElement *) (aux->FirstChild);
}
while (1);
if (j > n) // i.e. if (j >= n+1)
return -(n+1);
return i; // j == n, i.e. j is the index of the last leaf
}
else if (t == Element::OCTASECTED)
{
if (j == 0) return i;
return ((OctasectedElement *) elements[i])->Child[j-1];
}
}
if (j > 0)
return -1;
return i; // no refinement
}
void Mesh::BisectTriTrans(DenseMatrix &pointmat, Triangle *tri, int child)
{
double np[2];
if (child == 0) // left triangle
{
// Set the new coordinates of the vertices
np[0] = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
np[1] = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
pointmat(0,1) = pointmat(0,0); pointmat(1,1) = pointmat(1,0);
pointmat(0,0) = pointmat(0,2); pointmat(1,0) = pointmat(1,2);
pointmat(0,2) = np[0]; pointmat(1,2) = np[1];
}
else // right triangle
{
// Set the new coordinates of the vertices
np[0] = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
np[1] = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
pointmat(0,0) = pointmat(0,1); pointmat(1,0) = pointmat(1,1);
pointmat(0,1) = pointmat(0,2); pointmat(1,1) = pointmat(1,2);
pointmat(0,2) = np[0]; pointmat(1,2) = np[1];
}
}
void Mesh::BisectTetTrans(DenseMatrix &pointmat, Tetrahedron *tet, int child)
{
int i, j, redges[2], type, flag, ind[4];
double t[4];
tet->ParseRefinementFlag(redges, type, flag);
if (child == 0) // left tetrahedron
{
// Set the new coordinates of the vertices
pointmat(0,1) = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
pointmat(1,1) = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
pointmat(2,1) = 0.5 * ( pointmat(2,0) + pointmat(2,1) );
// Permute the vertices according to the type & redges
switch (redges[0])
{
case 2: ind[0] = 0; ind[1] = 3; ind[2] = 1; ind[3] = 2; break;
case 3: ind[0] = 1; ind[1] = 3; ind[2] = 2; ind[3] = 0; break;
case 5: ind[0] = 2; ind[1] = 3; ind[2] = 0; ind[3] = 1;
}
}
else // right tetrahedron
{
// Set the new coordinates of the vertices
pointmat(0,0) = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
pointmat(1,0) = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
pointmat(2,0) = 0.5 * ( pointmat(2,0) + pointmat(2,1) );
// Permute the vertices according to the type & redges
switch (redges[1])
{
case 1: ind[0] = 3; ind[1] = 1; ind[2] = 0; ind[3] = 2; break;
case 4: ind[0] = 3; ind[1] = 0; ind[2] = 2; ind[3] = 1; break;
case 5: ind[0] = 3; ind[1] = 2; ind[2] = 1; ind[3] = 0;
}
}
// Do the permutation
for (i = 0; i < 3; i++)
{
for (j = 0; j < 4; j++)
t[j] = pointmat(i,j);
for (j = 0; j < 4; j++)
pointmat(i,ind[j]) = t[j];
}
}
int Mesh::GetFineElemPath(int i, int j)
{
// if (Dim == 3)
{
if (elements[i]->GetType() == Element::BISECTED)
{
int n = 0, l = 0;
BisectedElement *aux = (BisectedElement *) elements[i];
do
{
int k = GetFineElemPath(aux->SecondChild, j-n);
if (k >= 0)
return ((k << 1)+1) << l;
n -= k; // (-k) is the number of the leaves in this SecondChild
// n is the number of the leaves in
// the SecondChild-ren so far
l++;
if (aux->FirstChild->GetType() != Element::BISECTED)
break;
aux = (BisectedElement *) (aux->FirstChild);
}
while (1);
if (j > n) // i.e. if (j >= n+1)
return -(n+1);
return 0; // j == n, i.e. j is the index of the last leaf
}
if (j > 0)
return -1;
}
return 0;
}
ElementTransformation * Mesh::GetFineElemTrans(int i, int j)
{
int t;
if (Dim == 2)
{
DenseMatrix &pm = Transformation.GetPointMat();
Transformation.Attribute = 0;
Transformation.ElementNo = 0;
t = elements[i]->GetType();
if (t == Element::QUADRISECTED)
{
t = ((QuadrisectedElement *)elements[i])->CoarseElem->GetType();
if (t == Element::QUADRILATERAL)
{
// quadrisected QUADRILATERAL
Transformation.SetFE(&QuadrilateralFE);
pm.SetSize(2, 4);
switch (j)
{
case 0:
pm(0,0) = 0.0; pm(1,0) = 0.0; // x; y;
pm(0,1) = 0.5; pm(1,1) = 0.0;
pm(0,2) = 0.5; pm(1,2) = 0.5;
pm(0,3) = 0.0; pm(1,3) = 0.5;
break;
case 1:
pm(0,0) = 0.5; pm(1,0) = 0.0;
pm(0,1) = 1.0; pm(1,1) = 0.0;
pm(0,2) = 1.0; pm(1,2) = 0.5;
pm(0,3) = 0.5; pm(1,3) = 0.5;
break;
case 2:
pm(0,0) = 0.5; pm(1,0) = 0.5;
pm(0,1) = 1.0; pm(1,1) = 0.5;
pm(0,2) = 1.0; pm(1,2) = 1.0;
pm(0,3) = 0.5; pm(1,3) = 1.0;
break;
case 3:
pm(0,0) = 0.0; pm(1,0) = 0.5;
pm(0,1) = 0.5; pm(1,1) = 0.5;
pm(0,2) = 0.5; pm(1,2) = 1.0;
pm(0,3) = 0.0; pm(1,3) = 1.0;
break;
default:
mfem_error("Mesh::GetFineElemTrans(...) 1");
}
}
else
{
// quadrisected TRIANGLE
Transformation.SetFE(&TriangleFE);
pm.SetSize(2, 3);
switch (j)
{
case 0:
pm(0,0) = 0.0; pm(0,1) = 0.5; pm(0,2) = 0.0; // x
pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 0.5; // y
break;
case 1:
pm(0,0) = 0.5; pm(0,1) = 1.0; pm(0,2) = 0.5;
pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 0.5;
break;
case 2:
pm(0,0) = 0.0; pm(0,1) = 0.5; pm(0,2) = 0.0;
pm(1,0) = 0.5; pm(1,1) = 0.5; pm(1,2) = 1.0;
break;
case 3:
pm(0,0) = 0.5; pm(0,1) = 0.0; pm(0,2) = 0.5;
pm(1,0) = 0.5; pm(1,1) = 0.5; pm(1,2) = 0.0;
break;
default:
mfem_error("Mesh::GetFineElemTrans(...) 2");
}
}
}
else if (t == Element::BISECTED)
{
// bisected TRIANGLE
Transformation.SetFE(&TriangleFE);
pm.SetSize(2, 3);
int path;
Element *E;
// pm is initialzed with the coordinates of the vertices of the
// reference triangle
pm(0,0) = 0.0; pm(0,1) = 1.0; pm(0,2) = 0.0;
pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 1.0;
path = GetFineElemPath(i, j);
E = elements[i];
while (E->GetType() == Element::BISECTED)
{
BisectedElement *aux = (BisectedElement *) E;
BisectTriTrans(pm, (Triangle *) aux->CoarseElem, path & 1);
E = (path & 1) ? elements[aux->SecondChild] : aux->FirstChild;
path = path >> 1;
}
}
else
{
// identity transformation
Transformation.SetFE(&TriangleFE);
pm.SetSize(2, 3);
pm(0,0) = 0.0; pm(0,1) = 1.0; pm(0,2) = 0.0;
pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 1.0;
}
return &Transformation;
}
else if (Dim == 3)
{
if (elements[i]->GetType() == Element::OCTASECTED)
{
int jj;
double dx, dy, dz;
DenseMatrix &pm = Transformation.GetPointMat();
Transformation.SetFE(&HexahedronFE);
Transformation.Attribute = 0;
Transformation.ElementNo = 0;
pm.SetSize(3, 8);
if (j < 4) dz = 0.0;
else dz = 0.5;
jj = j % 4;
if (jj < 2) dy = 0.0;
else dy = 0.5;
if (jj == 0 || jj == 3) dx = 0.0;
else dx = 0.5;
pm(0,0) = dx; pm(1,0) = dy; pm(2,0) = dz;
pm(0,1) = 0.5 + dx; pm(1,1) = dy; pm(2,1) = dz;
pm(0,2) = 0.5 + dx; pm(1,2) = 0.5 + dy; pm(2,2) = dz;
pm(0,3) = dx; pm(1,3) = 0.5 + dy; pm(2,3) = dz;
pm(0,4) = dx; pm(1,4) = dy; pm(2,4) = 0.5 + dz;
pm(0,5) = 0.5 + dx; pm(1,5) = dy; pm(2,5) = 0.5 + dz;
pm(0,6) = 0.5 + dx; pm(1,6) = 0.5 + dy; pm(2,6) = 0.5 + dz;
pm(0,7) = dx; pm(1,7) = 0.5 + dy; pm(2,7) = 0.5 + dz;
return &Transformation;
}
int path;
Element *E;
DenseMatrix &pm = Transformation.GetPointMat();
Transformation.SetFE(&TetrahedronFE);
Transformation.Attribute = 0;
Transformation.ElementNo = 0;
pm.SetSize(3, 4);
// pm is initialzed with the coordinates of the vertices of the
// reference tetrahedron
pm(0,0) = 0.0; pm(0,1) = 1.0; pm(0,2) = 0.0; pm(0,3) = 0.0;
pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 1.0; pm(1,3) = 0.0;
pm(2,0) = 0.0; pm(2,1) = 0.0; pm(2,2) = 0.0; pm(2,3) = 1.0;
path = GetFineElemPath(i, j);
E = elements[i];
while (E->GetType() == Element::BISECTED)
{
BisectedElement *aux = (BisectedElement *) E;
BisectTetTrans(pm, (Tetrahedron *) aux->CoarseElem, path & 1);
E = (path & 1) ? elements[aux->SecondChild] : aux->FirstChild;
path = path >> 1;
}
}
return &Transformation; // no refinement
}
void Mesh::PrintXG(ostream &out) const
{
int i, j;
Array<int> v;
if (Dim == 2)
{
// Print the type of the mesh.
if (Nodes == NULL)
out << "areamesh2\n\n";
else
out << "curved_areamesh2\n\n";
// Print the boundary elements.
out << NumOfBdrElements << '\n';
for (i = 0; i < NumOfBdrElements; i++)
{
boundary[i]->GetVertices(v);
out << boundary[i]->GetAttribute();
for (j = 0; j < v.Size(); j++)
out << ' ' << v[j] + 1;
out << '\n';
}
// Print the elements.
out << NumOfElements << '\n';
for (i = 0; i < NumOfElements; i++)
{
elements[i]->GetVertices(v);
out << elements[i]->GetAttribute() << ' ' << v.Size();
for (j = 0; j < v.Size(); j++)
out << ' ' << v[j] + 1;
out << '\n';
}
if (Nodes == NULL)
{
// Print the vertices.
out << NumOfVertices << '\n';
for (i = 0; i < NumOfVertices; i++)
{
out << vertices[i](0);
for (j = 1; j < Dim; j++)
out << ' ' << vertices[i](j);
out << '\n';
}
}
else
{
out << NumOfVertices << '\n';
Nodes->Save(out);
}
}
else // ===== Dim != 2 =====
{
if (Nodes)
{
mfem_error("Mesh::PrintXG(...) : Curved mesh in 3D");
}
if (meshgen == 1)
{
int nv;
const int *ind;
out << "NETGEN_Neutral_Format\n";
// print the vertices
out << NumOfVertices << '\n';
for (i = 0; i < NumOfVertices; i++)
{
for (j = 0; j < Dim; j++)
out << ' ' << vertices[i](j);
out << '\n';
}
// print the elements
out << NumOfElements << '\n';
for (i = 0; i < NumOfElements; i++)
{
nv = elements[i]->GetNVertices();
ind = elements[i]->GetVertices();
out << elements[i]->GetAttribute();
for (j = 0; j < nv; j++)
out << ' ' << ind[j]+1;
out << '\n';
}
// print the boundary information.
out << NumOfBdrElements << '\n';
for (i = 0; i < NumOfBdrElements; i++)
{
nv = boundary[i]->GetNVertices();
ind = boundary[i]->GetVertices();
out << boundary[i]->GetAttribute();
for (j = 0; j < nv; j++)
out << ' ' << ind[j]+1;
out << '\n';
}
}
else if (meshgen == 2) // TrueGrid
{
int nv;
const int *ind;
out << "TrueGrid\n"
<< "1 " << NumOfVertices << " " << NumOfElements
<< " 0 0 0 0 0 0 0\n"
<< "0 0 0 1 0 0 0 0 0 0 0\n"
<< "0 0 " << NumOfBdrElements << " 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
<< "0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
<< "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
for (i = 0; i < NumOfVertices; i++)
out << i+1 << " 0.0 " << vertices[i](0) << ' ' << vertices[i](1)
<< ' ' << vertices[i](2) << " 0.0\n";
for (i = 0; i < NumOfElements; i++)
{
nv = elements[i]->GetNVertices();
ind = elements[i]->GetVertices();
out << i+1 << ' ' << elements[i]->GetAttribute();
for (j = 0; j < nv; j++)
out << ' ' << ind[j]+1;
out << '\n';
}
for (i = 0; i < NumOfBdrElements; i++)
{
nv = boundary[i]->GetNVertices();
ind = boundary[i]->GetVertices();
out << boundary[i]->GetAttribute();
for (j = 0; j < nv; j++)
out << ' ' << ind[j]+1;
out << " 1.0 1.0 1.0 1.0\n";
}
}
}
out << flush;
}
void Mesh::Print(ostream &out) const
{
int i, j, nv, *v;
out << "MFEM mesh v1.0\n";
// optional
out <<
"\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
"# POINT = 0\n"
"# SEGMENT = 1\n"
"# TRIANGLE = 2\n"
"# SQUARE = 3\n"
"# TETRAHEDRON = 4\n"
"# CUBE = 5\n"
"#\n";
out << "\ndimension\n" << Dim
<< "\n\nelements\n" << NumOfElements << '\n';
for (i = 0; i < NumOfElements; i++)
{
out << elements[i]->GetAttribute() << ' '
<< elements[i]->GetGeometryType();
nv = elements[i]->GetNVertices();
v = elements[i]->GetVertices();
for (j = 0; j < nv; j++)
out << ' ' << v[j];
out << '\n';
}
out << "\nboundary\n" << NumOfBdrElements << '\n';
for (i = 0; i < NumOfBdrElements; i++)
{
out << boundary[i]->GetAttribute() << ' '
<< boundary[i]->GetGeometryType();
nv = boundary[i]->GetNVertices();
v = boundary[i]->GetVertices();
for (j = 0; j < nv; j++)
out << ' ' << v[j];
out << '\n';
}
out << "\nvertices\n" << NumOfVertices << '\n';
if (Nodes == NULL)
{
out << Dim << '\n';
for (i = 0; i < NumOfVertices; i++)
{
out << vertices[i](0);
for (j = 1; j < Dim; j++)
out << ' ' << vertices[i](j);
out << '\n';
}
}
else
{
out << "\nnodes\n";
Nodes->Save(out);
}
}
void Mesh::PrintVTK(ostream &out)
{
out <<
"# vtk DataFile Version 3.0\n"
"Generated by MFEM\n"
"ASCII\n"
"DATASET UNSTRUCTURED_GRID\n";
if (Nodes == NULL)
{
out << "POINTS " << NumOfVertices << " double\n";
for (int i = 0; i < NumOfVertices; i++)
{
out << vertices[i](0);
int j;
for (j = 1; j < Dim; j++)
out << ' ' << vertices[i](j);
for ( ; j < 3; j++)
out << ' ' << 0.0;
out << '\n';
}
}
else
{
Array<int> vdofs(3);
out << "POINTS " << Nodes->FESpace()->GetNDofs() << " double\n";
for (int i = 0; i < Nodes->FESpace()->GetNDofs(); i++)
{
vdofs.SetSize(1);
vdofs[0] = i;
Nodes->FESpace()->DofsToVDofs(vdofs);
out << (*Nodes)(vdofs[0]);
int j;
for (j = 1; j < Dim; j++)
out << ' ' << (*Nodes)(vdofs[j]);
for ( ; j < 3; j++)
out << ' ' << 0.0;
out << '\n';
}
}
int order = -1;
if (Nodes == NULL)
{
int size = 0;
for (int i = 0; i < NumOfElements; i++)
size += elements[i]->GetNVertices() + 1;
out << "CELLS " << NumOfElements << ' ' << size << '\n';
for (int i = 0; i < NumOfElements; i++)
{
const int *v = elements[i]->GetVertices();
const int nv = elements[i]->GetNVertices();
out << nv;
for (int j = 0; j < nv; j++)
out << ' ' << v[j];
out << '\n';
}
order = 1;
}
else
{
Array<int> dofs;
int size = 0;
for (int i = 0; i < NumOfElements; i++)
{
Nodes->FESpace()->GetElementDofs(i, dofs);
size += dofs.Size() + 1;
}
out << "CELLS " << NumOfElements << ' ' << size << '\n';
const char *fec_name = Nodes->FESpace()->FEColl()->Name();
if (!strcmp(fec_name, "Linear"))
order = 1;
else if (!strcmp(fec_name, "Quadratic"))
order = 2;
if (order == -1)
{
cerr << "Mesh::PrintVTK : can not save '"
<< fec_name << "' elements!" << endl;
mfem_error();
}
for (int i = 0; i < NumOfElements; i++)
{
Nodes->FESpace()->GetElementDofs(i, dofs);
out << dofs.Size();
if (order == 1)
{
for (int j = 0; j < dofs.Size(); j++)
out << ' ' << dofs[j];
}
else if (order == 2)
{
const int *vtk_mfem;
switch (elements[i]->GetGeometryType())
{
case Geometry::TRIANGLE:
case Geometry::SQUARE:
vtk_mfem = vtk_quadratic_hex; break; // identity map
case Geometry::TETRAHEDRON:
vtk_mfem = vtk_quadratic_tet; break;
case Geometry::CUBE:
vtk_mfem = vtk_quadratic_hex; break;
}
for (int j = 0; j < dofs.Size(); j++)
out << ' ' << dofs[vtk_mfem[j]];
}
out << '\n';
}
}
out << "CELL_TYPES " << NumOfElements << '\n';
for (int i = 0; i < NumOfElements; i++)
{
int vtk_cell_type;
if (order == 1)
{
switch (elements[i]->GetGeometryType())
{
case Geometry::TRIANGLE: vtk_cell_type = 5; break;
case Geometry::SQUARE: vtk_cell_type = 9; break;
case Geometry::TETRAHEDRON: vtk_cell_type = 10; break;
case Geometry::CUBE: vtk_cell_type = 12; break;
}
}
else if (order == 2)
{
switch (elements[i]->GetGeometryType())
{
case Geometry::TRIANGLE: vtk_cell_type = 22; break;
case Geometry::SQUARE: vtk_cell_type = 28; break;
case Geometry::TETRAHEDRON: vtk_cell_type = 24; break;
case Geometry::CUBE: vtk_cell_type = 29; break;
}
}
out << vtk_cell_type << '\n';
}
// write attributes
out << "CELL_DATA " << NumOfElements << '\n'
<< "SCALARS material int\n"
<< "LOOKUP_TABLE default\n";
for (int i = 0; i < NumOfElements; i++)
{
out << elements[i]->GetAttribute() << '\n';
}
}
void Mesh::PrintVTK(ostream &out, int ref)
{
int np, nc, size;
RefinedGeometry *RefG;
DenseMatrix pmat;
out <<
"# vtk DataFile Version 3.0\n"
"Generated by MFEM\n"
"ASCII\n"
"DATASET UNSTRUCTURED_GRID\n";
// count the points, cells, size
np = nc = size = 0;
for (int i = 0; i < GetNE(); i++)
{
int geom = GetElementBaseGeometry(i);
int nv = Geometries.GetVertices(geom)->GetNPoints();
RefG = GlobGeometryRefiner.Refine(geom, ref, 1);
np += RefG->RefPts.GetNPoints();
nc += RefG->RefGeoms.Size() / nv;
size += (RefG->RefGeoms.Size() / nv) * (nv + 1);
}
out << "POINTS " << np << " double\n";
// write the points
for (int i = 0; i < GetNE(); i++)
{
RefG = GlobGeometryRefiner.Refine(
GetElementBaseGeometry(i), ref, 1);
GetElementTransformation(i)->Transform(RefG->RefPts, pmat);
for (int j = 0; j < pmat.Width(); j++)
{
out << pmat(0, j) << ' ' << pmat(1, j) << ' ';
if (pmat.Height() == 2)
out << 0.0;
else
out << pmat(2, j);
out << '\n';
}
}
// write the cells
out << "CELLS " << nc << ' ' << size << '\n';
np = 0;
for (int i = 0; i < GetNE(); i++)
{
int geom = GetElementBaseGeometry(i);
int nv = Geometries.GetVertices(geom)->GetNPoints();
RefG = GlobGeometryRefiner.Refine(geom, ref, 1);
Array<int> &RG = RefG->RefGeoms;
for (int j = 0; j < RG.Size(); )
{
out << nv;
for (int k = 0; k < nv; k++, j++)
out << ' ' << np + RG[j];
out << '\n';
}
np += RefG->RefPts.GetNPoints();
}
out << "CELL_TYPES " << nc << '\n';
for (int i = 0; i < GetNE(); i++)
{
int geom = GetElementBaseGeometry(i);
int nv = Geometries.GetVertices(geom)->GetNPoints();
RefG = GlobGeometryRefiner.Refine(geom, ref, 1);
Array<int> &RG = RefG->RefGeoms;
int vtk_cell_type = 5;
switch (geom)
{
case Geometry::TRIANGLE: vtk_cell_type = 5; break;
case Geometry::SQUARE: vtk_cell_type = 9; break;
case Geometry::TETRAHEDRON: vtk_cell_type = 10; break;
case Geometry::CUBE: vtk_cell_type = 12; break;
}
for (int j = 0; j < RG.Size(); j += nv)
{
out << vtk_cell_type << '\n';
}
}
// write attributes (materials)
out << "CELL_DATA " << nc << '\n'
<< "SCALARS material int\n"
<< "LOOKUP_TABLE default\n";
for (int i = 0; i < GetNE(); i++)
{
int geom = GetElementBaseGeometry(i);
int nv = Geometries.GetVertices(geom)->GetNPoints();
RefG = GlobGeometryRefiner.Refine(geom, ref, 1);
int attr = GetAttribute(i);
for (int j = 0; j < RefG->RefGeoms.Size(); j += nv)
{
out << attr << '\n';
}
}
Array<int> coloring;
srandom(time(0));
double a = double(random()) / (double(RAND_MAX) + 1.);
int el0 = (int)floor(a * GetNE());
GetElementColoring(coloring, el0);
out << "SCALARS element_coloring int\n"
<< "LOOKUP_TABLE default\n";
for (int i = 0; i < GetNE(); i++)
{
int geom = GetElementBaseGeometry(i);
int nv = Geometries.GetVertices(geom)->GetNPoints();
RefG = GlobGeometryRefiner.Refine(geom, ref, 1);
for (int j = 0; j < RefG->RefGeoms.Size(); j += nv)
{
out << coloring[i] + 1 << '\n';
}
}
// prepare to write data
out << "POINT_DATA " << np << '\n';
}
void Mesh::GetElementColoring(Array<int> &colors, int el0)
{
int delete_el_to_el = (el_to_el) ? (0) : (1);
const Table &el_el = ElementToElementTable();
int num_el = GetNE(), stack_p, stack_top_p, max_num_col;
Array<int> el_stack(num_el);
const int *i_el_el = el_el.GetI();
const int *j_el_el = el_el.GetJ();
colors.SetSize(num_el);
colors = -2;
max_num_col = 1;
stack_p = stack_top_p = 0;
for (int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
{
if (colors[el] != -2)
continue;
colors[el] = -1;
el_stack[stack_top_p++] = el;
for ( ; stack_p < stack_top_p; stack_p++)
{
int i = el_stack[stack_p];
int num_nb = i_el_el[i+1] - i_el_el[i];
if (max_num_col < num_nb + 1)
max_num_col = num_nb + 1;
for (int j = i_el_el[i]; j < i_el_el[i+1]; j++)
{
int k = j_el_el[j];
if (colors[k] == -2)
{
colors[k] = -1;
el_stack[stack_top_p++] = k;
}
}
}
}
Array<int> col_marker(max_num_col);
for (stack_p = 0; stack_p < stack_top_p; stack_p++)
{
int i = el_stack[stack_p], col;
col_marker = 0;
for (int j = i_el_el[i]; j < i_el_el[i+1]; j++)
{
col = colors[j_el_el[j]];
if (col != -1)
col_marker[col] = 1;
}
for (col = 0; col < max_num_col; col++)
if (col_marker[col] == 0)
break;
colors[i] = col;
}
if (delete_el_to_el)
{
delete el_to_el;
el_to_el = NULL;
}
}
void Mesh::PrintWithPartitioning(int *partitioning, ostream &out) const
{
if (Dim != 3 && Dim != 2) return;
int i, j, k, l, nv, nbe, *v;
out << "MFEM mesh v1.0\n";
// optional
out <<
"\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
"# POINT = 0\n"
"# SEGMENT = 1\n"
"# TRIANGLE = 2\n"
"# SQUARE = 3\n"
"# TETRAHEDRON = 4\n"
"# CUBE = 5\n"
"#\n";
out << "\ndimension\n" << Dim
<< "\n\nelements\n" << NumOfElements << '\n';
for (i = 0; i < NumOfElements; i++)
{
out << elements[i]->GetAttribute() << ' '
<< elements[i]->GetGeometryType();
nv = elements[i]->GetNVertices();
v = elements[i]->GetVertices();
for (j = 0; j < nv; j++)
out << ' ' << v[j];
out << '\n';
}
nbe = 0;
for (i = 0; i < NumOfFaces; i++)
{
if ((l = faces_info[i].Elem2No) >= 0)
{
k = partitioning[faces_info[i].Elem1No];
l = partitioning[l];
if (k != l)
nbe += 2;
}
else
nbe++;
}
out << "\nboundary\n" << nbe << '\n';
for (i = 0; i < NumOfFaces; i++)
{
if ((l = faces_info[i].Elem2No) >= 0)
{
k = partitioning[faces_info[i].Elem1No];
l = partitioning[l];
if (k != l)
{
nv = faces[i]->GetNVertices();
v = faces[i]->GetVertices();
out << k+1 << ' ' << faces[i]->GetGeometryType();
for (j = 0; j < nv; j++)
out << ' ' << v[j];
out << '\n';
out << l+1 << ' ' << faces[i]->GetGeometryType();
for (j = nv-1; j >= 0; j--)
out << ' ' << v[j];
out << '\n';
}
}
else
{
k = partitioning[faces_info[i].Elem1No];
nv = faces[i]->GetNVertices();
v = faces[i]->GetVertices();
out << k+1 << ' ' << faces[i]->GetGeometryType();
for (j = 0; j < nv; j++)
out << ' ' << v[j];
out << '\n';
}
}
out << "\nvertices\n" << NumOfVertices << '\n';
if (Nodes == NULL)
{
out << Dim << '\n';
for (i = 0; i < NumOfVertices; i++)
{
out << vertices[i](0);
for (j = 1; j < Dim; j++)
out << ' ' << vertices[i](j);
out << '\n';
}
}
else
{
out << "\nnodes\n";
Nodes->Save(out);
}
}
void Mesh::PrintElementsWithPartitioning(int *partitioning,
ostream &out,
int interior_faces)
{
if (Dim != 3 && Dim != 2) return;
int i, j, k, l, s;
int nv;
const int *ind;
int *vcount = new int[NumOfVertices];
for (i = 0; i < NumOfVertices; i++)
vcount[i] = 0;
for (i = 0; i < NumOfElements; i++)
{
nv = elements[i]->GetNVertices();
ind = elements[i]->GetVertices();
for (j = 0; j < nv; j++)
vcount[ind[j]]++;
}
int *voff = new int[NumOfVertices+1];
voff[0] = 0;
for (i = 1; i <= NumOfVertices; i++)
voff[i] = vcount[i-1] + voff[i-1];
int **vown = new int*[NumOfVertices];
for (i = 0; i < NumOfVertices; i++)
vown[i] = new int[vcount[i]];
// 2D
if (Dim == 2)
{
int nv, nbe;
int *ind;
Table edge_el;
Transpose(ElementToEdgeTable(), edge_el);
// Fake printing of the elements.
for (i = 0; i < NumOfElements; i++)
{
nv = elements[i]->GetNVertices();
ind = elements[i]->GetVertices();
for (j = 0; j < nv; j++)
{
vcount[ind[j]]--;
vown[ind[j]][vcount[ind[j]]] = i;
}
}
for (i = 0; i < NumOfVertices; i++)
vcount[i] = voff[i+1] - voff[i];
nbe = 0;
for (i = 0; i < edge_el.Size(); i++)
{
const int *el = edge_el.GetRow(i);
if (edge_el.RowSize(i) > 1)
{
k = partitioning[el[0]];
l = partitioning[el[1]];
if (interior_faces || k != l)
nbe += 2;
}
else
nbe++;
}
// Print the type of the mesh and the boundary elements.
out << "areamesh2\n\n" << nbe << '\n';
for (i = 0; i < edge_el.Size(); i++)
{
const int *el = edge_el.GetRow(i);
if (edge_el.RowSize(i) > 1)
{
k = partitioning[el[0]];
l = partitioning[el[1]];
if (interior_faces || k != l)
{
Array<int> ev;
GetEdgeVertices(i,ev);
out << k+1; // attribute
for (j = 0; j < 2; j++)
for (s = 0; s < vcount[ev[j]]; s++)
if (vown[ev[j]][s] == el[0])
out << ' ' << voff[ev[j]]+s+1;
out << '\n';
out << l+1; // attribute
for (j = 1; j >= 0; j--)
for (s = 0; s < vcount[ev[j]]; s++)
if (vown[ev[j]][s] == el[1])
out << ' ' << voff[ev[j]]+s+1;
out << '\n';
}
}
else
{
k = partitioning[el[0]];
Array<int> ev;
GetEdgeVertices(i,ev);
out << k+1; // attribute
for (j = 0; j < 2; j++)
for (s = 0; s < vcount[ev[j]]; s++)
if (vown[ev[j]][s] == el[0])
out << ' ' << voff[ev[j]]+s+1;
out << '\n';
}
}
// Print the elements.
out << NumOfElements << '\n';
for (i = 0; i < NumOfElements; i++)
{
nv = elements[i]->GetNVertices();
ind = elements[i]->GetVertices();
out << partitioning[i]+1 << ' '; // use subdomain number as attribute
out << nv << ' ';
for (j = 0; j < nv; j++)
{
out << ' ' << voff[ind[j]]+vcount[ind[j]]--;
vown[ind[j]][vcount[ind[j]]] = i;
}
out << '\n';
}
for (i = 0; i < NumOfVertices; i++)
vcount[i] = voff[i+1] - voff[i];
// Print the vertices.
out << voff[NumOfVertices] << '\n';
for (i = 0; i < NumOfVertices; i++)
for (k = 0; k < vcount[i]; k++)
{
for (j = 0; j < Dim; j++)
out << vertices[i](j) << ' ';
out << '\n';
}
out << flush;
return;
}
// Dim is 3
if (meshgen == 1)
{
out << "NETGEN_Neutral_Format\n";
// print the vertices
out << voff[NumOfVertices] << '\n';
for (i = 0; i < NumOfVertices; i++)
for (k = 0; k < vcount[i]; k++)
{
for (j = 0; j < Dim; j++)
out << ' ' << vertices[i](j);
out << '\n';
}
// print the elements
out << NumOfElements << '\n';
for (i = 0; i < NumOfElements; i++)
{
nv = elements[i]->GetNVertices();
ind = elements[i]->GetVertices();
out << partitioning[i]+1; // use subdomain number as attribute
for (j = 0; j < nv; j++)
{
out << ' ' << voff[ind[j]]+vcount[ind[j]]--;
vown[ind[j]][vcount[ind[j]]] = i;
}
out << '\n';
}
for (i = 0; i < NumOfVertices; i++)
vcount[i] = voff[i+1] - voff[i];
// print the boundary information.
int k, l, nbe;
nbe = 0;
for (i = 0; i < NumOfFaces; i++)
if ((l = faces_info[i].Elem2No) >= 0)
{
k = partitioning[faces_info[i].Elem1No];
l = partitioning[l];
if (interior_faces || k != l)
nbe += 2;
}
else
nbe++;
out << nbe << '\n';
for (i = 0; i < NumOfFaces; i++)
if ((l = faces_info[i].Elem2No) >= 0)
{
k = partitioning[faces_info[i].Elem1No];
l = partitioning[l];
if (interior_faces || k != l)
{
nv = faces[i]->GetNVertices();
ind = faces[i]->GetVertices();
out << k+1; // attribute
for (j = 0; j < nv; j++)
for (s = 0; s < vcount[ind[j]]; s++)
if (vown[ind[j]][s] == faces_info[i].Elem1No)
out << ' ' << voff[ind[j]]+s+1;
out << '\n';
out << l+1; // attribute
for (j = nv-1; j >= 0; j--)
for (s = 0; s < vcount[ind[j]]; s++)
if (vown[ind[j]][s] == faces_info[i].Elem2No)
out << ' ' << voff[ind[j]]+s+1;
out << '\n';
}
}
else
{
k = partitioning[faces_info[i].Elem1No];
nv = faces[i]->GetNVertices();
ind = faces[i]->GetVertices();
out << k+1; // attribute
for (j = 0; j < nv; j++)
for (s = 0; s < vcount[ind[j]]; s++)
if (vown[ind[j]][s] == faces_info[i].Elem1No)
out << ' ' << voff[ind[j]]+s+1;
out << '\n';
}
for (i = 0; i < NumOfVertices; i++)
delete [] vown[i];
}
else if (meshgen == 2) // TrueGrid
{
// count the number of the boundary elements.
int k, l, nbe;
nbe = 0;
for (i = 0; i < NumOfFaces; i++)
if ((l = faces_info[i].Elem2No) >= 0)
{
k = partitioning[faces_info[i].Elem1No];
l = partitioning[l];
if (interior_faces || k != l)
nbe += 2;
}
else
nbe++;
out << "TrueGrid\n"
<< "1 " << voff[NumOfVertices] << " " << NumOfElements
<< " 0 0 0 0 0 0 0\n"
<< "0 0 0 1 0 0 0 0 0 0 0\n"
<< "0 0 " << nbe << " 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
<< "0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
<< "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
for (i = 0; i < NumOfVertices; i++)
for (k = 0; k < vcount[i]; k++)
out << voff[i]+k << " 0.0 " << vertices[i](0) << ' '
<< vertices[i](1) << ' ' << vertices[i](2) << " 0.0\n";
for (i = 0; i < NumOfElements; i++)
{
nv = elements[i]->GetNVertices();
ind = elements[i]->GetVertices();
out << i+1 << ' ' << partitioning[i]+1; // partitioning as attribute
for (j = 0; j < nv; j++)
{
out << ' ' << voff[ind[j]]+vcount[ind[j]]--;
vown[ind[j]][vcount[ind[j]]] = i;
}
out << '\n';
}
for (i = 0; i < NumOfVertices; i++)
vcount[i] = voff[i+1] - voff[i];
// boundary elements
for (i = 0; i < NumOfFaces; i++)
if ((l = faces_info[i].Elem2No) >= 0)
{
k = partitioning[faces_info[i].Elem1No];
l = partitioning[l];
if (interior_faces || k != l)
{
nv = faces[i]->GetNVertices();
ind = faces[i]->GetVertices();
out << k+1; // attribute
for (j = 0; j < nv; j++)
for (s = 0; s < vcount[ind[j]]; s++)
if (vown[ind[j]][s] == faces_info[i].Elem1No)
out << ' ' << voff[ind[j]]+s+1;
out << " 1.0 1.0 1.0 1.0\n";
out << l+1; // attribute
for (j = nv-1; j >= 0; j--)
for (s = 0; s < vcount[ind[j]]; s++)
if (vown[ind[j]][s] == faces_info[i].Elem2No)
out << ' ' << voff[ind[j]]+s+1;
out << " 1.0 1.0 1.0 1.0\n";
}
}
else
{
k = partitioning[faces_info[i].Elem1No];
nv = faces[i]->GetNVertices();
ind = faces[i]->GetVertices();
out << k+1; // attribute
for (j = 0; j < nv; j++)
for (s = 0; s < vcount[ind[j]]; s++)
if (vown[ind[j]][s] == faces_info[i].Elem1No)
out << ' ' << voff[ind[j]]+s+1;
out << " 1.0 1.0 1.0 1.0\n";
}
}
out << flush;
delete [] vcount;
delete [] voff;
delete [] vown;
}
void Mesh::ScaleSubdomains(double sf)
{
int i,j,k;
Array<int> vert;
DenseMatrix pointmat;
int na = attributes.Size();
double *cg = new double[na*Dim];
int *nbea = new int[na];
int *vn = new int[NumOfVertices];
for (i = 0; i < NumOfVertices; i++)
vn[i] = 0;
for (i = 0; i < na; i++)
{
for (j = 0; j < Dim; j++)
cg[i*Dim+j] = 0.0;
nbea[i] = 0;
}
for (i = 0; i < NumOfElements; i++)
{
GetElementVertices(i, vert);
for (k = 0; k < vert.Size(); k++)
vn[vert[k]] = 1;
}
for (i = 0; i < NumOfElements; i++)
{
int bea = GetAttribute(i)-1;
GetPointMatrix(i, pointmat);
GetElementVertices(i, vert);
for (k = 0; k < vert.Size(); k++)
if (vn[vert[k]] == 1)
{
nbea[bea]++;
for (j = 0; j < Dim; j++)
cg[bea*Dim+j] += pointmat(j,k);
vn[vert[k]] = 2;
}
}
for (i = 0; i < NumOfElements; i++)
{
int bea = GetAttribute(i)-1;
GetElementVertices (i, vert);
for (k = 0; k < vert.Size(); k++)
if (vn[vert[k]])
{
for (j = 0; j < Dim; j++)
vertices[vert[k]](j) = sf*vertices[vert[k]](j) +
(1-sf)*cg[bea*Dim+j]/nbea[bea];
vn[vert[k]] = 0;
}
}
delete [] cg;
delete [] nbea;
delete [] vn;
}
void Mesh::ScaleElements(double sf)
{
int i,j,k;
Array<int> vert;
DenseMatrix pointmat;
int na = NumOfElements;
double *cg = new double[na*Dim];
int *nbea = new int[na];
int *vn = new int[NumOfVertices];
for (i = 0; i < NumOfVertices; i++)
vn[i] = 0;
for (i = 0; i < na; i++)
{
for (j = 0; j < Dim; j++)
cg[i*Dim+j] = 0.0;
nbea[i] = 0;
}
for (i = 0; i < NumOfElements; i++)
{
GetElementVertices(i, vert);
for (k = 0; k < vert.Size(); k++)
vn[vert[k]] = 1;
}
for (i = 0; i < NumOfElements; i++)
{
int bea = i;
GetPointMatrix(i, pointmat);
GetElementVertices(i, vert);
for (k = 0; k < vert.Size(); k++)
if (vn[vert[k]] == 1)
{
nbea[bea]++;
for (j = 0; j < Dim; j++)
cg[bea*Dim+j] += pointmat(j,k);
vn[vert[k]] = 2;
}
}
for (i = 0; i < NumOfElements; i++)
{
int bea = i;
GetElementVertices(i, vert);
for (k = 0; k < vert.Size(); k++)
if (vn[vert[k]])
{
for (j = 0; j < Dim; j++)
vertices[vert[k]](j) = sf*vertices[vert[k]](j) +
(1-sf)*cg[bea*Dim+j]/nbea[bea];
vn[vert[k]] = 0;
}
}
delete [] cg;
delete [] nbea;
delete [] vn;
}
void Mesh::Transform(void (*f)(const Vector&, Vector&))
{
if (Nodes == NULL)
{
Vector vold(Dim), vnew(NULL, Dim);
for (int i = 0; i < vertices.Size(); i++)
{
for (int j = 0; j < Dim; j++)
vold(j) = vertices[i](j);
vnew.SetData(vertices[i]());
(*f)(vold, vnew);
}
}
else
{
GridFunction xnew(Nodes->FESpace());
VectorFunctionCoefficient f_pert(Dim, f);
xnew.ProjectCoefficient(f_pert);
*Nodes = xnew;
}
}
void Mesh::FreeElement(Element *E)
{
#ifdef MFEM_USE_MEMALLOC
if (E)
switch (E->GetType())
{
case Element::TETRAHEDRON: TetMemory.Free((Tetrahedron *)E); break;
case Element::BISECTED: BEMemory.Free((BisectedElement *)E); break;
default: delete E; break;
}
#else
delete E;
#endif
}
Mesh::~Mesh()
{
int i;
if (own_nodes) delete Nodes;
for (i = 0; i < NumOfElements; i++)
FreeElement(elements[i]);
for (i = 0; i < NumOfBdrElements; i++)
FreeElement(boundary[i]);
for (i = 0; i < faces.Size(); i++)
FreeElement(faces[i]);
DeleteTables();
}
Reference in New Issue View Git Blame Copy Permalink