Select Git revision
-
Christophe Geuzaine authored
tutorial/ -> tutorials/
Christophe Geuzaine authoredtutorial/ -> tutorials/
GModelVertexArrays.cpp 15.41 KiB
// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
#include <cmath>
#include "GmshMessage.h"
#include "GmshDefines.h"
#include "GModel.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
#include "MTrihedron.h"
#include "MElementCut.h"
#include "Context.h"
#include "VertexArray.h"
#include "SmoothData.h"
static const double curvedRepTol = 1.e-5;
unsigned int getColorByEntity(GEntity *e)
{
if(e->getSelection()) { // selection
return CTX::instance()->color.geom.selection;
}
else if(e->useColor()) { // forced from a script
return e->getColor();
}
else if(CTX::instance()->mesh.colorCarousel == 1) { // by elementary entity
return CTX::instance()->color.mesh.carousel[abs(e->tag() % 20)];
}
else if(CTX::instance()->mesh.colorCarousel == 2) { // by physical entity
int np = e->physicals.size();
int p = np ? e->physicals[np - 1] : 0;
return CTX::instance()->color.mesh.carousel[abs(p % 20)];
}
else {
return CTX::instance()->color.fg;
}
}
static unsigned int getColorByElement(MElement *ele)
{
if(ele->getVisibility() > 1) { // selection
return CTX::instance()->color.geom.selection;
}
else if(CTX::instance()->mesh.colorCarousel == 0) { // by element type
switch(ele->getType()) {
case TYPE_LIN: return CTX::instance()->color.mesh.line;
case TYPE_TRI: return CTX::instance()->color.mesh.triangle;
case TYPE_QUA: return CTX::instance()->color.mesh.quadrangle;
case TYPE_TET: return CTX::instance()->color.mesh.tetrahedron;
case TYPE_HEX: return CTX::instance()->color.mesh.hexahedron;
case TYPE_PRI: return CTX::instance()->color.mesh.prism;
case TYPE_PYR: return CTX::instance()->color.mesh.pyramid;
case TYPE_TRIH: return CTX::instance()->color.mesh.trihedron;
default: return CTX::instance()->color.mesh.node;
}
}
else if(CTX::instance()->mesh.colorCarousel == 3) { // by partition
return CTX::instance()
->color.mesh.carousel[std::abs(ele->getPartition() % 20)];
}
else {
// by elementary or physical entity (this is not perfect (since
// e.g. a triangle can have no vertices categorized on a surface),
// but it's the best we can do "fast" since we don't store the
// associated entity in the element
for(std::size_t i = 0; i < ele->getNumVertices(); i++) {
GEntity *e = ele->getVertex(i)->onWhat();
if(e && (e->dim() == ele->getDim())) return getColorByEntity(e);
}
}
return CTX::instance()->color.fg;
}
static double evalClipPlane(int clip, double x, double y, double z)
{
return CTX::instance()->clipPlane[clip][0] * x +
CTX::instance()->clipPlane[clip][1] * y +
CTX::instance()->clipPlane[clip][2] * z +
CTX::instance()->clipPlane[clip][3];
}
static double intersectClipPlane(int clip, MElement *ele)
{
MVertex *v = ele->getVertex(0);
double val = evalClipPlane(clip, v->x(), v->y(), v->z());
for(std::size_t i = 1; i < ele->getNumVertices(); i++) {
v = ele->getVertex(i);
if(val * evalClipPlane(clip, v->x(), v->y(), v->z()) <= 0)
return 0.; // the element intersects the cut plane
}
return val;
}
bool isElementVisible(MElement *ele)
{
if(!ele->getVisibility()) return false;
if(CTX::instance()->mesh.qualitySup) {
double q;
if(CTX::instance()->mesh.qualityType == 3)
q = ele->distoShapeMeasure();
else if(CTX::instance()->mesh.qualityType == 2)
q = ele->gammaShapeMeasure();
else if(CTX::instance()->mesh.qualityType == 1)
q = ele->minSIGEShapeMeasure();
else
q = ele->minSICNShapeMeasure();
if(q < CTX::instance()->mesh.qualityInf ||
q > CTX::instance()->mesh.qualitySup)
return false;
}
if(CTX::instance()->mesh.radiusSup) {
double r = ele->maxEdge();
if(r < CTX::instance()->mesh.radiusInf ||
r > CTX::instance()->mesh.radiusSup)
return false;
}
if(CTX::instance()->clipWholeElements) {
bool hidden = false;
for(int clip = 0; clip < 6; clip++) {
if(CTX::instance()->mesh.clip & (1 << clip)) {
if(ele->getDim() < 3 && CTX::instance()->clipOnlyVolume) {
}
else {
double d = intersectClipPlane(clip, ele);
if(ele->getDim() == 3 &&
CTX::instance()->clipOnlyDrawIntersectingVolume && d) {
hidden = true;
break;
}
else if(d < 0) {
hidden = true;
break;
}
}
}
}
if(hidden) return false;
}
return true;
}
template <class T> static bool areAllElementsVisible(std::vector<T *> &elements)
{
for(unsigned int i = 0; i < elements.size(); i++)
if(!isElementVisible(elements[i])) return false;
return true;
}
template <class T> static bool areSomeElementsCurved(std::vector<T *> &elements)
{
for(unsigned int i = 0; i < elements.size(); i++)
if(elements[i]->getPolynomialOrder() > 1) return true;
return false;
}
template <class T>
static void addSmoothNormals(GEntity *e, std::vector<T *> &elements)
{
for(unsigned int i = 0; i < elements.size(); i++) {
MElement *ele = elements[i];
const bool curved =
(ele->getPolynomialOrder() > 1) &&
(ele->maxDistToStraight() > curvedRepTol * ele->getInnerRadius());
SPoint3 pc(0., 0., 0.);
if(CTX::instance()->mesh.explode != 1.) pc = ele->barycenter();
for(int j = 0; j < ele->getNumFacesRep(curved); j++) {
double x[3], y[3], z[3];
SVector3 n[3];
ele->getFaceRep(curved, j, x, y, z, n);
for(int k = 0; k < 3; k++) {
if(CTX::instance()->mesh.explode != 1.) {
x[k] = pc[0] + CTX::instance()->mesh.explode * (x[k] - pc[0]);
y[k] = pc[1] + CTX::instance()->mesh.explode * (y[k] - pc[1]);
z[k] = pc[2] + CTX::instance()->mesh.explode * (z[k] - pc[2]);
}
e->model()->normals->add(x[k], y[k], z[k], n[k][0], n[k][1], n[k][2]);
}
}
}
}
template <class T>
static void addElementsInArrays(GEntity *e, std::vector<T *> &elements,
bool edges, bool faces)
{
#if defined(_OPENMP)
#pragma omp parallel for schedule(dynamic)
#endif
for(unsigned int i = 0; i < elements.size(); i++) {
MElement *ele = elements[i];
if(!isElementVisible(ele) || ele->getDim() < 1) continue;
unsigned int c = getColorByElement(ele);
unsigned int col[4] = {c, c, c, c};
const bool curved =
(ele->getPolynomialOrder() > 1) &&
(ele->maxDistToStraight() > curvedRepTol * ele->getInnerRadius());
SPoint3 pc(0., 0., 0.);
if(CTX::instance()->mesh.explode != 1.) pc = ele->barycenter();
if(edges) {
bool unique = e->dim() > 1 && !CTX::instance()->pickElements;
for(int j = 0; j < ele->getNumEdgesRep(curved); j++) {
double x[2], y[2], z[2];
SVector3 n[2];
ele->getEdgeRep(curved, j, x, y, z, n);
if(CTX::instance()->mesh.explode != 1.) {
for(int k = 0; k < 2; k++) {
x[k] = pc[0] + CTX::instance()->mesh.explode * (x[k] - pc[0]);
y[k] = pc[1] + CTX::instance()->mesh.explode * (y[k] - pc[1]);
z[k] = pc[2] + CTX::instance()->mesh.explode * (z[k] - pc[2]);
}
}
if(e->dim() == 2 && CTX::instance()->mesh.smoothNormals)
for(int k = 0; k < 2; k++)
e->model()->normals->get(x[k], y[k], z[k], n[k][0], n[k][1],
n[k][2]);
#if defined(_OPENMP)
#pragma omp critical
#endif
{
e->va_lines->add(x, y, z, n, col, ele, unique);
}
}
}
if(faces) {
bool unique = e->dim() > 2 && !CTX::instance()->pickElements;
bool skin = e->dim() > 2 && CTX::instance()->mesh.drawSkinOnly;
for(int j = 0; j < ele->getNumFacesRep(curved); j++) {
double x[3], y[3], z[3];
SVector3 n[3];
ele->getFaceRep(curved, j, x, y, z, n);
if(CTX::instance()->mesh.explode != 1.) {
for(int k = 0; k < 3; k++) {
x[k] = pc[0] + CTX::instance()->mesh.explode * (x[k] - pc[0]);
y[k] = pc[1] + CTX::instance()->mesh.explode * (y[k] - pc[1]);
z[k] = pc[2] + CTX::instance()->mesh.explode * (z[k] - pc[2]);
}
}
if(e->dim() == 2 && CTX::instance()->mesh.smoothNormals)
for(int k = 0; k < 3; k++)
e->model()->normals->get(x[k], y[k], z[k], n[k][0], n[k][1],
n[k][2]);
#if defined(_OPENMP)
#pragma omp critical
#endif
{
e->va_triangles->add(x, y, z, n, col, ele, unique, skin);
}
}
}
}
}
class initMeshGEdge {
private:
int _estimateNumLines(GEdge *e)
{
int num = 0;
if(CTX::instance()->mesh.lines) {
num += e->lines.size();
if(areSomeElementsCurved(e->lines)) num *= 2;
}
return num + 100;
}
public:
void operator()(GEdge *e)
{
e->deleteVertexArrays();
if(!e->getVisibility()) return;
e->setAllElementsVisible(CTX::instance()->mesh.lines &&
areAllElementsVisible(e->lines));
if(CTX::instance()->mesh.lines) {
e->va_lines = new VertexArray(2, _estimateNumLines(e));
addElementsInArrays(e, e->lines, CTX::instance()->mesh.lines, false);
e->va_lines->finalize();
}
}
};
class initSmoothNormalsGFace {
public:
void operator()(GFace *f)
{
addSmoothNormals(f, f->triangles);
addSmoothNormals(f, f->quadrangles);
addSmoothNormals(f, f->polygons);
}
};
class initMeshGFace {
private:
bool _curved;
int _estimateNumLines(GFace *f)
{
int num = 0;
if(CTX::instance()->mesh.surfacesEdges) {
num += (3 * f->triangles.size() + 4 * f->quadrangles.size() +
4 * f->polygons.size()) /
2;
if(CTX::instance()->mesh.explode != 1.) num *= 2;
if(_curved) num *= 2;
}
return num + 100;
}
int _estimateNumTriangles(GFace *f)
{
int num = 0;
if(CTX::instance()->mesh.surfacesFaces) {
num += (f->triangles.size() + 2 * f->quadrangles.size() +
2 * f->polygons.size());
if(_curved) num *= 4;
}
return num + 100;
}
public:
void operator()(GFace *f)
{
f->deleteVertexArrays();
if(!f->getVisibility()) return;
f->setAllElementsVisible(CTX::instance()->mesh.triangles &&
areAllElementsVisible(f->triangles) &&
CTX::instance()->mesh.quadrangles &&
areAllElementsVisible(f->quadrangles));
bool edg = CTX::instance()->mesh.surfacesEdges;
bool fac = CTX::instance()->mesh.surfacesFaces;
if(edg || fac) {
_curved = (areSomeElementsCurved(f->triangles) ||
areSomeElementsCurved(f->quadrangles));
f->va_lines = new VertexArray(2, edg ? _estimateNumLines(f) : 100);
f->va_triangles = new VertexArray(3, fac ? _estimateNumTriangles(f) : 100);
if(CTX::instance()->mesh.triangles)
addElementsInArrays(f, f->triangles, edg, fac);
if(CTX::instance()->mesh.quadrangles)
addElementsInArrays(f, f->quadrangles, edg, fac);
addElementsInArrays(f, f->polygons, edg, fac);
f->va_lines->finalize();
f->va_triangles->finalize();
}
}
};
class initMeshGRegion {
private:
bool _curved;
int _estimateIfClipped(int num)
{
if(CTX::instance()->clipWholeElements) {
for(int clip = 0; clip < 6; clip++) {
if(CTX::instance()->mesh.clip & (1 << clip)){
if(CTX::instance()->clipOnlyDrawIntersectingVolume){
// let be more aggressive than num^{2/3}
return (int)sqrt((double)num);
}
else{
// why not :-)
return num / 4;
}
}
}
}
return num;
}
int _estimateNumLines(GRegion *r)
{
int num = 0;
if(CTX::instance()->mesh.volumesEdges) {
// suppose edge shared by 4 elements on averge (pessmistic)
int numLP = 0;
for(unsigned int i = 0; i < r->polyhedra.size(); i++)
numLP += 2 * r->polyhedra[i]->getNumEdges();
num += (12 * r->tetrahedra.size() + 24 * r->hexahedra.size() +
18 * r->prisms.size() + 16 * r->pyramids.size() +
10 * r->trihedra.size() + numLP) /
4;
num = _estimateIfClipped(num);
if(CTX::instance()->mesh.explode != 1.) num *= 4;
if(_curved) num *= 2;
}
return num + 100;
}
int _estimateNumTriangles(GRegion *r)
{
int num = 0;
if(CTX::instance()->mesh.volumesFaces) {
int numFP = 0;
for(unsigned int i = 0; i < r->polyhedra.size(); i++)
numFP += r->polyhedra[i]->getNumFaces();
num += (4 * r->tetrahedra.size() + 12 * r->hexahedra.size() +
8 * r->prisms.size() + 6 * r->pyramids.size() +
4 * r->trihedra.size() + numFP) /
2;
num = _estimateIfClipped(num);
if(CTX::instance()->mesh.explode != 1.) num *= 2;
if(_curved) num *= 4;
}
return num + 100;
}
public:
void operator()(GRegion *r)
{
r->deleteVertexArrays();
if(!r->getVisibility()) return;
r->setAllElementsVisible(
CTX::instance()->mesh.tetrahedra &&
areAllElementsVisible(r->tetrahedra) && CTX::instance()->mesh.hexahedra &&
areAllElementsVisible(r->hexahedra) && CTX::instance()->mesh.prisms &&
areAllElementsVisible(r->prisms) && CTX::instance()->mesh.pyramids &&
areAllElementsVisible(r->pyramids) && CTX::instance()->mesh.trihedra &&
areAllElementsVisible(r->trihedra));
bool edg = CTX::instance()->mesh.volumesEdges;
bool fac = CTX::instance()->mesh.volumesFaces;
if(edg || fac) {
_curved = (areSomeElementsCurved(r->tetrahedra) ||
areSomeElementsCurved(r->hexahedra) ||
areSomeElementsCurved(r->prisms) ||
areSomeElementsCurved(r->pyramids) ||
areSomeElementsCurved(r->trihedra));
r->va_lines = new VertexArray(2, edg ? _estimateNumLines(r) : 100);
r->va_triangles = new VertexArray(3, fac ? _estimateNumTriangles(r) : 100);
if(CTX::instance()->mesh.tetrahedra)
addElementsInArrays(r, r->tetrahedra, edg, fac);
if(CTX::instance()->mesh.hexahedra)
addElementsInArrays(r, r->hexahedra, edg, fac);
if(CTX::instance()->mesh.prisms)
addElementsInArrays(r, r->prisms, edg, fac);
if(CTX::instance()->mesh.pyramids)
addElementsInArrays(r, r->pyramids, edg, fac);
if(CTX::instance()->mesh.trihedra)
addElementsInArrays(r, r->trihedra, edg, fac);
addElementsInArrays(r, r->polyhedra, edg, fac);
r->va_lines->finalize();
r->va_triangles->finalize();
}
}
};
bool GModel::fillVertexArrays()
{
if(!getVisibility() || !CTX::instance()->mesh.changed) return false;
Msg::Debug("Mesh has changed: reinitializing vertex arrays");
int status = getMeshStatus();
if(status >= 1 && CTX::instance()->mesh.changed & ENT_CURVE)
std::for_each(firstEdge(), lastEdge(), initMeshGEdge());
if(status >= 2 && CTX::instance()->mesh.changed & ENT_SURFACE) {
if(normals) delete normals;
normals = new smooth_normals(CTX::instance()->mesh.angleSmoothNormals);
if(CTX::instance()->mesh.smoothNormals)
std::for_each(firstFace(), lastFace(), initSmoothNormalsGFace());
std::for_each(firstFace(), lastFace(), initMeshGFace());
}
if(status >= 3 && CTX::instance()->mesh.changed & ENT_VOLUME)
std::for_each(firstRegion(), lastRegion(), initMeshGRegion());
return true;
}