implement meshing from 2D geometry

fix vtk writer for 2d meshes
extend vtk writer to write also domain index per cell
This commit is contained in:
Bryn Lloyd 2019-03-01 22:21:55 +01:00
parent 353795062c
commit b19a17f2d7
3 changed files with 112 additions and 38 deletions

View File

@ -50,64 +50,68 @@ void WriteVtkFormat (const Mesh & mesh,
outfile << "POINTS " << np << " double\n";
for (int i=0; i<np; i++)
{
auto & p = mesh.Point(i+1);
auto & p = mesh.Point(i+1);
outfile << p[0] << " " << p[1] << " " << p[2] << "\n";
}
std::vector<int> types;
std::vector<int> domains;
if (ne > 0)
{
unsigned int size = 0;
for (int i=0; i<ne; i++)
size += mesh.VolumeElement(i+1).GetNV() + 1; // only save "linear" corners
size += mesh.VolumeElement(i+1).GetNV() + 1; // only save "linear" corners
outfile << "CELLS " << ne << " " << size << "\n";
for (int i=0; i<ne; i++)
{
auto& el = mesh.VolumeElement(i+1);
switch (el.GetType())
{
case TET:
case TET10: // reorder to follow VTK convention & zero based indices
auto& el = mesh.VolumeElement(i+1);
domains.push_back(el.GetIndex());
switch (el.GetType())
{
case TET:
case TET10: // reorder to follow VTK convention & zero based indices
outfile << 4 << " " << el[0]-1 << " " << el[1]-1 << " " << el[3]-1 << " " << el[2]-1 << "\n";
types.push_back(10);
break;
case PRISM: // reorder to follow VTK convention & zero based indices
types.push_back(10);
break;
case PRISM: // reorder to follow VTK convention & zero based indices
outfile << 6 << " "
<< el[0]-1 << " " << el[2]-1 << " " << el[1]-1 << " "
<< el[3]-1 << " " << el[5]-1 << " " << el[4]-1 << "\n";
types.push_back(13);
break;
default:
throw ngcore::Exception("Unexpected element type");
break;
}
<< el[0]-1 << " " << el[2]-1 << " " << el[1]-1 << " "
<< el[3]-1 << " " << el[5]-1 << " " << el[4]-1 << "\n";
types.push_back(13);
break;
default:
throw ngcore::Exception("Unexpected element type");
break;
}
}
}
else
{
unsigned int size = 0;
for (int i=0; i<ne; i++)
size += mesh.SurfaceElement(i+1).GetNV() + 1;
for (int i=0; i<nse; i++)
size += mesh.SurfaceElement(i+1).GetNV() + 1;
outfile << "CELLS " << nse << " " << size << "\n";
for (int i=0; i<nse; i++)
{
auto& el = mesh.SurfaceElement(i+1);
switch (el.GetType())
{
case TRIG:
case TRIG6:
outfile << el[0]-1 << " " << el[1]-1 << " " << el[2]-1 << "\n";
types.push_back(5);
break;
case QUAD:
outfile << el[0]-1 << " " << el[1]-1 << " " << el[2]-1 << " " << el[3]-1 << "\n";
types.push_back(9);
default:
throw ngcore::Exception("Unexpected element type");
break;
}
auto& el = mesh.SurfaceElement(i+1);
domains.push_back(el.GetIndex());
switch (el.GetType())
{
case TRIG:
case TRIG6:
outfile << 3 << " " << el[0]-1 << " " << el[1]-1 << " " << el[2]-1 << "\n";
types.push_back(5);
break;
case QUAD:
outfile << 4 << " " << el[0]-1 << " " << el[1]-1 << " " << el[2]-1 << " " << el[3]-1 << "\n";
types.push_back(9);
break;
default:
throw ngcore::Exception("Unexpected element type");
break;
}
}
}
@ -117,6 +121,14 @@ void WriteVtkFormat (const Mesh & mesh,
outfile << type_id << "\n";
}
outfile << "CELL_DATA " << domains.size() << "\n";
outfile << "SCALARS scalars int 1\n";
outfile << "LOOKUP_TABLE default\n";
for (auto id: domains)
{
outfile << id << "\n";
}
outfile.close();
}
}

View File

@ -591,6 +591,68 @@ namespace nglib
}
DLL_HEADER Ng_Geometry_2D * Ng_NewGeometry_2D ()
{
SplineGeometry2d * geom = new SplineGeometry2d();
return (Ng_Geometry_2D *)geom;
}
DLL_HEADER void Ng_DeleteGeometry_2D (Ng_Geometry_2D * geom)
{
if (geom)
{
SplineGeometry2d* spline_geom = (SplineGeometry2d*)geom;
delete spline_geom;
geom = NULL;
}
}
DLL_HEADER void Ng_AppendPoint_2D (Ng_Geometry_2D* geom, double * x, double h)
{
if (geom)
{
SplineGeometry2d* spline_geom = (SplineGeometry2d*)geom;
Point<2> p(x[0],x[1]);
spline_geom->AppendPoint(p, h);
}
}
DLL_HEADER void Ng_AppendLineSegment_2D (Ng_Geometry_2D* geom, int n1, int n2,
int leftdomain, int rightdomain, double h)
{
if (geom)
{
SplineGeometry2d* spline_geom = (SplineGeometry2d*)geom;
// zero-offset!
LineSeg<2>* line = new LineSeg<2>(spline_geom->geompoints[n1-1], spline_geom->geompoints[n2-1]);
SplineSegExt* seg = new SplineSegExt(*line);
seg->leftdom = leftdomain;
seg->rightdom = rightdomain;
seg->hmax = h;
seg->copyfrom = -1;
seg->bc = 1;
spline_geom->AppendSegment(seg);
}
}
DLL_HEADER void Ng_AppendSplinSegment_2D (Ng_Geometry_2D* geom, int n1, int n2, int n3,
int leftdomain, int rightdomain, double h)
{
if (geom)
{
SplineGeometry2d* spline_geom = (SplineGeometry2d*)geom;
// zero-offset!
SplineSeg3<2>* line = new SplineSeg3<2>(spline_geom->geompoints[n1-1], spline_geom->geompoints[n2-1], spline_geom->geompoints[n3-1]);
SplineSegExt* seg = new SplineSegExt(*line);
seg->leftdom = leftdomain;
seg->rightdom = rightdomain;
seg->hmax = h;
seg->copyfrom = -1;
seg->bc = 1;
spline_geom->AppendSegment(seg);
}
}
DLL_HEADER Ng_Result Ng_GenerateMesh_2D (Ng_Geometry_2D * geom,
Ng_Mesh ** mesh,
Ng_Meshing_Parameters * mp)
@ -601,9 +663,6 @@ namespace nglib
shared_ptr<Mesh> m(new Mesh, &NOOP_Deleter);
MeshFromSpline2D (*(SplineGeometry2d*)geom, m, mparam);
// new shared_ptr<Mesh> (m); // hack to keep mesh m alive
cout << m->GetNSE() << " elements, " << m->GetNP() << " points" << endl;
*mesh = (Ng_Mesh*)m.get();
return NG_OK;

View File

@ -625,6 +625,9 @@ DLL_HEADER void Ng_AppendPoint_2D(Ng_Geometry_2D* geom, double * x, double h);
DLL_HEADER void Ng_AppendLineSegment_2D(Ng_Geometry_2D* geom, int pi1, int pi2,
int leftdomain, int rightdomain, double h);
DLL_HEADER void Ng_AppendSplinSegment_2D(Ng_Geometry_2D* geom, int pi1, int pi2, int pi3,
int leftdomain, int rightdomain, double h);
DLL_HEADER void Ng_AddPoint_2D (Ng_Mesh * mesh, double * x);
DLL_HEADER void Ng_AddBoundarySeg_2D (Ng_Mesh * mesh, int pi1, int pi2, int domain_in = -1, int domain_out = -1);