Select Git revision
GRegion.cpp 13.91 KiB
// Gmsh - Copyright (C) 1997-2016 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.
#include <sstream>
#include "GModel.h"
#include "GRegion.h"
#include "GFace.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
#include "MTrihedron.h"
#include "MElementCut.h"
#include "GmshMessage.h"
#include "VertexArray.h"
#include "boundaryLayersData.h"
#include "GmshDefines.h"
GRegion::GRegion(GModel *model, int tag)
: GEntity (model, tag), compound(0)
{
resetMeshAttributes();
}
GRegion::~GRegion()
{
std::list<GFace*>::iterator it = l_faces.begin();
while(it != l_faces.end()){
(*it)->delRegion(this);
++it;
}
deleteMesh();
}
void GRegion::deleteMesh()
{
for(unsigned int i = 0; i < mesh_vertices.size(); i++) delete mesh_vertices[i];
mesh_vertices.clear();
transfinite_vertices.clear();
for(unsigned int i = 0; i < tetrahedra.size(); i++) delete tetrahedra[i];
tetrahedra.clear();
for(unsigned int i = 0; i < hexahedra.size(); i++) delete hexahedra[i];
hexahedra.clear();
for(unsigned int i = 0; i < prisms.size(); i++) delete prisms[i];
prisms.clear();
for(unsigned int i = 0; i < pyramids.size(); i++) delete pyramids[i];
pyramids.clear();
for(unsigned int i = 0; i < trihedra.size(); i++) delete trihedra[i];
trihedra.clear();
for(unsigned int i = 0; i < polyhedra.size(); i++) delete polyhedra[i];
polyhedra.clear();
deleteVertexArrays();
model()->destroyMeshCaches();
}
unsigned int GRegion::getNumMeshElements()
{
return tetrahedra.size() + hexahedra.size() + prisms.size() + pyramids.size() +
trihedra.size() + polyhedra.size();
}
unsigned int GRegion::getNumMeshParentElements()
{
unsigned int n = 0;
for(unsigned int i = 0; i < polyhedra.size(); i++)
if(polyhedra[i]->ownsParent())
n++;
return n;
}
void GRegion::getNumMeshElements(unsigned *const c) const
{
c[0] += tetrahedra.size();
c[1] += hexahedra.size();
c[2] += prisms.size();
c[3] += pyramids.size();
c[4] += trihedra.size();
c[5] += polyhedra.size();
}
MElement *const *GRegion::getStartElementType(int type) const
{
switch(type) {
case 0:
if(tetrahedra.empty()) return 0; // msvc would throw an exception
return reinterpret_cast<MElement *const *>(&tetrahedra[0]);
case 1:
if(hexahedra.empty()) return 0; // msvc would throw an exception
return reinterpret_cast<MElement *const *>(&hexahedra[0]);
case 2:
if(prisms.empty()) return 0; // msvc would throw an exception
return reinterpret_cast<MElement *const *>(&prisms[0]);
case 3:
if(pyramids.empty()) return 0; // msvc would throw an exception
return reinterpret_cast<MElement *const *>(&pyramids[0]);
case 4:
if(trihedra.empty()) return 0;
return reinterpret_cast<MElement *const *>(&trihedra[0]);
case 5:
if(polyhedra.empty()) return 0;
return reinterpret_cast<MElement *const *>(&polyhedra[0]);
}
return 0;
}
MElement *GRegion::getMeshElement(unsigned int index)
{
if(index < tetrahedra.size())
return tetrahedra[index];
else if(index < tetrahedra.size() + hexahedra.size())
return hexahedra[index - tetrahedra.size()];
else if(index < tetrahedra.size() + hexahedra.size() + prisms.size())
return prisms[index - tetrahedra.size() - hexahedra.size()];
else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
pyramids.size())
return pyramids[index - tetrahedra.size() - hexahedra.size() - prisms.size()];
else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
pyramids.size() + trihedra.size())
return trihedra[index - tetrahedra.size() - hexahedra.size() - prisms.size() -
pyramids.size()];
else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
pyramids.size() + trihedra.size() + polyhedra.size())
return polyhedra[index - tetrahedra.size() - hexahedra.size() - prisms.size() -
pyramids.size() - trihedra.size()];
return 0;
}
void GRegion::resetMeshAttributes()
{
meshAttributes.recombine3D = 0;
meshAttributes.method = MESH_UNSTRUCTURED;
meshAttributes.extrude = 0;
meshAttributes.QuadTri = NO_QUADTRI;
}
SBoundingBox3d GRegion::bounds() const
{
SBoundingBox3d res;
if(geomType() != DiscreteVolume){
std::list<GFace*>::const_iterator it = l_faces.begin();
for(; it != l_faces.end(); it++)
res += (*it)->bounds();
}
else{
for(unsigned int i = 0; i < mesh_vertices.size(); i++)
res += mesh_vertices[i]->point();
}
return res;
}
SOrientedBoundingBox GRegion::getOBB()
{
if (!_obb) {
std::vector<SPoint3> vertices;
std::list<GFace*> b_faces = faces();
for (std::list<GFace*>::iterator b_face = b_faces.begin();
b_face != b_faces.end(); b_face++) {
if((*b_face)->getNumMeshVertices() > 0) {
int N = (*b_face)->getNumMeshVertices();
for (int i = 0; i < N; i++) {
MVertex* mv = (*b_face)->getMeshVertex(i);
vertices.push_back(mv->point());
}
std::list<GEdge*> eds = (*b_face)->edges();
for(std::list<GEdge*>::iterator ed = eds.begin(); ed != eds.end(); ed++) {
int N2 = (*ed)->getNumMeshVertices();
for (int i = 0; i < N2; i++) {
MVertex* mv = (*ed)->getMeshVertex(i);
vertices.push_back(mv->point());
}
// Don't forget to add the first and last vertices...
SPoint3 pt1((*ed)->getBeginVertex()->x(),
(*ed)->getBeginVertex()->y(),
(*ed)->getBeginVertex()->z());
SPoint3 pt2((*ed)->getEndVertex()->x(),
(*ed)->getEndVertex()->y(),
(*ed)->getEndVertex()->z());
vertices.push_back(pt1);
vertices.push_back(pt2);
}
}
else if ((*b_face)->buildSTLTriangulation()) {
for (unsigned int i = 0; i < (*b_face)->stl_vertices.size(); i++){
GPoint p = (*b_face)->point((*b_face)->stl_vertices[i]);
vertices.push_back(SPoint3(p.x(), p.y(), p.z()));
}
}
else {
int N = 10;
std::list<GEdge*> b_edges = (*b_face)->edges();
for (std::list<GEdge*>::iterator b_edge = b_edges.begin();
b_edge != b_edges.end(); b_edge++) {
Range<double> tr = (*b_edge)->parBounds(0);
for (int j = 0; j < N; j++) {
double t = tr.low() + (double)j / (double)(N - 1) * (tr.high() - tr.low());
GPoint p = (*b_edge)->point(t);
SPoint3 pt(p.x(), p.y(), p.z());
vertices.push_back(pt);
}
}
}
}
_obb = SOrientedBoundingBox::buildOBB(vertices);
}
return SOrientedBoundingBox(_obb);
}
void GRegion::setVisibility(char val, bool recursive)
{
GEntity::setVisibility(val);
if(recursive){
std::list<GFace*>::iterator it = l_faces.begin();
while(it != l_faces.end()){
(*it)->setVisibility(val, recursive);
++it;
}
}
}
std::string GRegion::getAdditionalInfoString()
{
std::ostringstream sstream;
if(l_faces.size()){
sstream << "{";
for(std::list<GFace*>::iterator it = l_faces.begin(); it != l_faces.end(); ++it){
if(it != l_faces.begin()) sstream << " ";
sstream << (*it)->tag();
}
sstream << "}";
}
if(meshAttributes.method == MESH_TRANSFINITE)
sstream << " transfinite";
if(meshAttributes.extrude)
sstream << " extruded";
return sstream.str();
}
void GRegion::writeGEO(FILE *fp)
{
if(geomType() == DiscreteVolume) return;
if(l_faces.size()){
fprintf(fp, "Surface Loop(%d) = ", tag());
for(std::list<GFace*>::iterator it = l_faces.begin(); it != l_faces.end(); it++) {
if(it != l_faces.begin())
fprintf(fp, ", %d", (*it)->tag());
else
fprintf(fp, "{%d", (*it)->tag());
}
fprintf(fp, "};\n");
fprintf(fp, "Volume(%d) = {%d};\n", tag(), tag());
}
if(meshAttributes.method == MESH_TRANSFINITE){
fprintf(fp, "Transfinite Volume {%d}", tag());
if(meshAttributes.corners.size()){
fprintf(fp, " = {");
for(unsigned int i = 0; i < meshAttributes.corners.size(); i++){
if(i) fprintf(fp, ",");
fprintf(fp, "%d", meshAttributes.corners[i]->tag());
}
fprintf(fp, "}");
}
fprintf(fp, ";\n");
if(meshAttributes.QuadTri != NO_QUADTRI )
fprintf(fp, "TransfQuadTri {%d};\n", tag());
}
}
std::list<GEdge*> GRegion::edges() const
{
std::list<GEdge*> e;
std::list<GFace*>::const_iterator it = l_faces.begin();
while(it != l_faces.end()){
std::list<GEdge*> e2;
e2 = (*it)->edges();
std::list<GEdge*>::const_iterator it2 = e2.begin();
while (it2 != e2.end()){
GEdge *edge = *it2;
// FIXME: we need to fix the compound design and decide what to do; same
// thing for faces() (either store or compute the entities, either use
// original or compound entities, etc.)
if(edge->getCompound())
edge = (GEdge*)edge->getCompound();
if(std::find(e.begin(), e.end(), edge) == e.end())
e.push_back(edge);
++it2;
}
++it;
}
return e;
}
bool GRegion::edgeConnected(GRegion *r) const
{
std::list<GEdge*> e = edges();
std::list<GEdge*> e2 = r->edges();
std::list<GEdge*>::const_iterator it = e.begin();
while(it != e.end()){
if(std::find(e2.begin(), e2.end(), *it) != e2.end())
return true;
++it;
}
return false;
}
void GRegion::replaceFaces(std::list<GFace*> &new_faces)
{
replaceFacesInternal(new_faces);
if (l_faces.size() != new_faces.size()){
Msg::Error("Impossible to replace faces in region %d (%d vs %d)",
tag(), l_faces.size(), new_faces.size());
}
std::list<GFace*>::iterator it = l_faces.begin();
std::list<GFace*>::iterator it2 = new_faces.begin();
std::list<int>::iterator it3 = l_dirs.begin();
std::list<int> newdirs;
for ( ; it != l_faces.end(); ++it, ++it2, ++it3){
(*it)->delRegion(this);
(*it2)->addRegion(this);
// if ((*it2)->getBeginVertex() == (*it)->getBeginVertex())
newdirs.push_back(*it3);
// else
// newdirs.push_back(-(*it3));
}
l_faces = new_faces;
l_dirs = newdirs;
}
double GRegion::computeSolidProperties(std::vector<double> cg,
std::vector<double> inertia)
{
std::list<GFace*>::iterator it = l_faces.begin();
std::list<int>::iterator itdir = l_dirs.begin();
double volumex = 0;
double volumey = 0;
double volumez = 0;
double surface = 0;
cg[0] = cg[1] = cg[2] = 0.0;
for ( ; it != l_faces.end(); ++it,++itdir){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *e = (*it)->triangles[i];
int npt;
IntPt *pts;
e->getIntegrationPoints (2*(e->getPolynomialOrder()-1)+3, &npt, &pts);
for (int j=0;j<npt;j++){
SPoint3 pt;
// compute x,y,z of the integration point
e->pnt(pts[j].pt[0], pts[j].pt[1], pts[j].pt[2], pt);
double jac[3][3];
// compute normal
double detJ = e->getJacobian(pts[j].pt[0], pts[j].pt[1], pts[j].pt[2], jac);
SVector3 n(jac[2][0], jac[2][1], jac[2][2]);
n.normalize();
n *= (double)*itdir;
surface += detJ* pts[j].weight;
volumex += detJ * n.x() * pt.x() * pts[j].weight;
volumey += detJ * n.y() * pt.y() * pts[j].weight;
volumez += detJ * n.z() * pt.z() * pts[j].weight;
cg[0] += detJ * n.x() * (pt.x() * pt.x()) * pts[j].weight * 0.5;
cg[1] += detJ * n.y() * (pt.y() * pt.y()) * pts[j].weight * 0.5;
cg[2] += detJ * n.z() * (pt.z() * pt.z()) * pts[j].weight * 0.5;
}
}
}
printf("%g -- %g %g %g\n", surface, volumex, volumey, volumez);
double volume = volumex;
cg[0] /= volume;
cg[1] /= volume;
cg[2] /= volume;
it = l_faces.begin();
itdir = l_dirs.begin();
inertia[0] = inertia[1] = inertia[2] = inertia[3] = inertia[4] = inertia[5] = 0.0;
for ( ; it != l_faces.end(); ++it,++itdir){
for (unsigned int i = 0; i < (*it)->getNumMeshElements(); ++i){
MElement *e = (*it)->getMeshElement(i);
int npt;
IntPt *pts;
e->getIntegrationPoints(2 * (e->getPolynomialOrder() - 1) + 3, &npt, &pts);
for (int j = 0; j < npt; j++){
SPoint3 pt;
// compute x,y,z of the integration point
e->pnt(pts[j].pt[0], pts[j].pt[1], pts[j].pt[2], pt);
double jac[3][3];
// compute normal
double detJ = e->getJacobian(pts[j].pt[0], pts[j].pt[1], pts[j].pt[2], jac);
SVector3 n(jac[2][0], jac[2][1], jac[2][2]);
n *= (double)*itdir;
inertia[0] += pts[j].weight * detJ * n.x() *
(pt.x() - cg[0]) * (pt.x() - cg[0]) * (pt.x() - cg[0]) / 3.0;
inertia[1] += pts[j].weight * detJ * n.y() *
(pt.y() - cg[1]) * (pt.y() - cg[1]) * (pt.y() - cg[1]) / 3.0;
inertia[2] += pts[j].weight * detJ * n.z() *
(pt.z() - cg[2]) * (pt.z() - cg[2]) * (pt.z() - cg[2]) / 3.0;
inertia[3] += pts[j].weight * detJ * n.x() *
(pt.y() - cg[1]) * (pt.x() - cg[0]) * (pt.x() - cg[0]) / 3.0;
inertia[4] += pts[j].weight * detJ * n.x() *
(pt.z() - cg[2]) * (pt.x() - cg[0]) * (pt.x() - cg[0]) / 3.0;
inertia[5] += pts[j].weight * detJ * n.y() *
(pt.z() - cg[2]) * (pt.y() - cg[1]) * (pt.y() - cg[1]) / 3.0;
}
}
}
return volume;
}
std::list<GVertex*> GRegion::vertices() const
{
std::set<GVertex*> v;
for (std::list<GFace*>::const_iterator it = l_faces.begin(); it != l_faces.end() ; ++it){
const GFace *gf = *it;
std::list<GVertex*> vs = gf->vertices();
v.insert(vs.begin(), vs.end());
}
std::list<GVertex*> res;
res.insert(res.begin(), v.begin(), v.end());
return res;
}
void GRegion::addElement(int type, MElement *e)
{
switch (type){
case TYPE_TET:
addTetrahedron((MTetrahedron*) e);
break;
case TYPE_HEX:
addHexahedron((MHexahedron*) e);
break;
case TYPE_PRI:
addPrism((MPrism*) e);
break;
case TYPE_PYR:
addPyramid((MPyramid*) e);
break;
case TYPE_TRIH:
addTrihedron((MTrihedron*) e);
break;
case TYPE_POLYH:
addPolyhedron((MPolyhedron*) e);
break;
default:
Msg::Error("Trying to add unsupported element in region");
}
}