From ad2b36e6e7f5bf66d7fe128d5a60dd50bf1614ba Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-Fran=C3=A7ois=20Remacle=20=28students=29?=
<jean-francois.remacle@uclouvain.be>
Date: Mon, 3 Aug 2009 14:42:30 +0000
Subject: [PATCH] fix display problems (gbricteux)
---
Geo/MElementCut.cpp | 96 ++---
Graphics/Makefile | 1 +
Graphics/drawMesh.cpp | 29 +-
contrib/DiscreteIntegration/DILevelset.cpp | 390 +++++++++---------
contrib/DiscreteIntegration/DILevelset.h | 95 +++--
contrib/DiscreteIntegration/Integration3D.cpp | 25 +-
6 files changed, 331 insertions(+), 305 deletions(-)
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index d7919e6f72..588a3b1a27 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -290,12 +290,19 @@ void assignPhysicals(GModel *GM, std::vector<int> gePhysicals, int reg, int dim,
physicals[dim][reg][gePhysicals[i]] = GM->getPhysicalName(dim, gePhysicals[i]);
}
+bool equalV(MVertex *v, DI_Point p)
+{
+ return (fabs(v->x() - p.x()) < 1.e-15 && fabs(v->y() - p.y()) < 1.e-15 &&
+ fabs(v->z() - p.z()) < 1.e-15);
+}
+
static void elementCutMesh (MElement *e, gLevelset *ls, GEntity *ge, GModel *GM, int &numEle,
std::map<int, MVertex*> &vertexMap,
std::vector<MVertex*> &newVertices,
std::map<int, std::vector<MElement*> > elements[10],
std::map<int, std::vector<MElement*> > border[2],
- std::map<int, std::map<int, std::string> > physicals[4])
+ std::map<int, std::map<int, std::string> > physicals[4],
+ std::map<int, std::vector<GEntity*> > &entityCut)
{
int elementary = ge->tag();
int eType = e->getTypeForMSH();
@@ -403,10 +410,7 @@ static void elementCutMesh (MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
if(numV == -1) {
unsigned int k;
for(k = 0; k < newVertices.size(); k++)
- if(subTetras[i].x(j) == newVertices[k]->x() &&
- subTetras[i].y(j) == newVertices[k]->y() &&
- subTetras[i].z(j) == newVertices[k]->z())
- break;
+ if(equalV(newVertices[k], subTetras[i].pt(j))) break;
if(k == newVertices.size()) {
mv[j] = new MVertex(subTetras[i].x(j), subTetras[i].y(j),
subTetras[i].z(j), 0, numV);
@@ -444,10 +448,7 @@ static void elementCutMesh (MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
if(numV == -1) {
unsigned int k;
for(k = 0; k < newVertices.size(); k++)
- if(surfTriangles[i].x(j) == newVertices[k]->x() &&
- surfTriangles[i].y(j) == newVertices[k]->y() &&
- surfTriangles[i].z(j) == newVertices[k]->z())
- break;
+ if(equalV(newVertices[k], surfTriangles[i].pt(j))) break;
if(k == newVertices.size()) {
mv[j] = new MVertex(surfTriangles[i].x(j), surfTriangles[i].y(j),
surfTriangles[i].z(j), 0, numV);
@@ -468,7 +469,7 @@ static void elementCutMesh (MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
MTriangleBorder *tri = new MTriangleBorder(mv[0], mv[1], mv[2],
p1, p2, ++numEle, ePart);
border[1][surfTriangles[i].lsTag()].push_back(tri);
- assignPhysicals(GM, gePhysicals, surfTriangles[i].lsTag(), 2, physicals);
+ entityCut[surfTriangles[i].lsTag()].push_back(ge);
}
}
@@ -527,10 +528,7 @@ static void elementCutMesh (MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
if(numV == -1) {
unsigned int k;
for(k = 0; k < newVertices.size(); k++)
- if(subTriangles[i].x(j) == newVertices[k]->x() &&
- subTriangles[i].y(j) == newVertices[k]->y() &&
- subTriangles[i].z(j) == newVertices[k]->z())
- break;
+ if(equalV(newVertices[k], subTriangles[i].pt(j))) break;
if(k == newVertices.size()) {
mv[j] = new MVertex(subTriangles[i].x(j), subTriangles[i].y(j),
subTriangles[i].z(j), 0, numV);
@@ -568,10 +566,7 @@ static void elementCutMesh (MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
if(numV == -1) {
unsigned int k;
for(k = 0; k < newVertices.size(); k++)
- if(boundLines[i].x(j) == newVertices[k]->x() &&
- boundLines[i].y(j) == newVertices[k]->y() &&
- boundLines[i].z(j) == newVertices[k]->z())
- break;
+ if(equalV(newVertices[k], boundLines[i].pt(j))) break;
if(k == newVertices.size()) {
mv[j] = new MVertex(boundLines[i].x(j), boundLines[i].y(j),
boundLines[i].z(j), 0, numV);
@@ -591,7 +586,7 @@ static void elementCutMesh (MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
}
MLineBorder *lin = new MLineBorder(mv[0], mv[1], p1, p2, ++numEle, ePart);
border[0][boundLines[i].lsTag()].push_back(lin);
- assignPhysicals(GM, gePhysicals, boundLines[i].lsTag(), 1, physicals);
+ entityCut[boundLines[i].lsTag()].push_back(ge);
}
}
@@ -623,10 +618,7 @@ static void elementCutMesh (MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
if(numV == -1) {
unsigned int k;
for(k = 0; k < newVertices.size(); k++)
- if(lines[i].x(j) == newVertices[k]->x() &&
- lines[i].y(j) == newVertices[k]->y() &&
- lines[i].z(j) == newVertices[k]->z())
- break;
+ if(equalV(newVertices[k], lines[i].pt(j))) break;
if(k == newVertices.size()) {
mv[j] = new MVertex(lines[i].x(j), lines[i].y(j), lines[i].z(j), 0, numV);
newVertices.push_back(mv[j]);
@@ -683,8 +675,9 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
std::map<int, std::vector<MElement*> > border[2];
std::vector<MVertex*> newVertices;
-
std::vector<GEntity*> gmEntities;
+ std::map<int, std::vector<GEntity*> > entityCut;
+
gm->getEntities(gmEntities);
int numEle = gm->getNumMeshElements();
@@ -692,11 +685,42 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
for(int j = 0; j < gmEntities[i]->getNumMeshElements(); j++) {
MElement *e = gmEntities[i]->getMeshElement(j);
elementCutMesh (e, ls, gmEntities[i], gm, numEle,
- vertexMap, newVertices, elements, border, physicals);
+ vertexMap, newVertices, elements, border, physicals, entityCut);
cutGM->getMeshPartitions().insert(e->getPartition());
}
}
+ // add borders in elements and change the tag if it's already used
+ std::map<int, std::vector<MElement*> >::iterator itbo, itel;
+ for(itbo = border[0].begin(); itbo != border[0].end(); itbo++) {
+ int reg = itbo->first;
+ if(elements[1].count(reg)) {
+ itel = elements[1].end(); itel--;
+ reg = itel->first + 1;
+ }
+ for(unsigned int j = 0; j < itbo->second.size(); j++)
+ elements[1][reg].push_back(itbo->second[j]);
+ std::map<int, std::vector<GEntity*> >::iterator itge = entityCut.find(itbo->first);
+ for(unsigned int j = 0; j < itge->second.size(); j++)
+ assignPhysicals(gm, itge->second[j]->physicals, reg, 1, physicals);
+ }
+ for(itbo = border[1].begin(); itbo != border[1].end(); itbo++) {
+ int reg = itbo->first;
+ if(elements[2].count(reg)) {
+ itel = elements[2].end(); itel--;
+ reg = itel->first + 1;
+ }
+ if(elements[3].count(reg)) {
+ itel = elements[3].end(); itel--;
+ reg = std::max(reg, itel->first + 1);
+ }
+ for(unsigned int j = 0; j < itbo->second.size(); j++)
+ elements[2][reg].push_back(itbo->second[j]);
+ std::map<int, std::vector<GEntity*> >::iterator itge = entityCut.find(itbo->first);
+ for(unsigned int j = 0; j < itge->second.size(); j++)
+ assignPhysicals(gm, itge->second[j]->physicals, reg, 2, physicals);
+ }
+
// number the new vertices and add in vertexMap
std::map<int, MVertex* >::iterator itend = vertexMap.end(); itend--;
int num = itend->first;
@@ -705,28 +729,6 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
vertexMap[num] = newVertices[i];
}printf("numbering vertices finished : %d vertices \n",vertexMap.size());
- // add borders in elements and check if the tag is not used
- std::map<int, std::vector<MElement*> >::iterator it = border[0].begin();
- for(; it != border[0].end(); it++) {
- int n = it->first;
- if(elements[1].find(n) != elements[1].end())
- n = elements[1].end()->first;
- for(int j = 0; j < it->second.size(); j++)
- elements[1][n].push_back(it->second[j]);
- it->second.clear();
- }
- it = border[1].begin();
- for(; it != border[1].end(); it++) {
- int n = it->first;
- if(elements[2].find(n) != elements[2].end())
- n = elements[2].end()->first;
- if(elements[3].find(n) != elements[3].end())
- n = std::max(n, elements[3].end()->first);
- for(int j = 0; j < it->second.size(); j++)
- elements[2][n].push_back(it->second[j]);
- it->second.clear();
- }
-
return cutGM;
#else
Msg::Error("Gmsh need to be compiled with levelset support to cut mesh");
diff --git a/Graphics/Makefile b/Graphics/Makefile
index c4853528e4..8942f888f9 100644
--- a/Graphics/Makefile
+++ b/Graphics/Makefile
@@ -98,6 +98,7 @@ drawMesh${OBJEXT}: drawMesh.cpp drawContext.h ../Geo/SBoundingBox3d.h \
../Geo/MElement.h ../Geo/MTetrahedron.h ../Geo/MElement.h \
../Geo/MHexahedron.h ../Geo/MElement.h ../Geo/MPrism.h \
../Geo/MElement.h ../Geo/MPyramid.h ../Geo/MElement.h \
+ ../Geo/MElementCut.h ../Geo/MElement.h \
../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
../Common/OS.h gl2ps.h ../Common/VertexArray.h ../Common/SmoothData.h \
../Post/PView.h ../Post/PViewData.h
diff --git a/Graphics/drawMesh.cpp b/Graphics/drawMesh.cpp
index 2127c97554..08a2487ad7 100644
--- a/Graphics/drawMesh.cpp
+++ b/Graphics/drawMesh.cpp
@@ -15,6 +15,7 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
+#include "MElementCut.h"
#include "Context.h"
#include "OS.h"
#include "gl2ps.h"
@@ -212,7 +213,7 @@ static void drawNormals(drawContext *ctx, std::vector<T*> &elements)
{
glColor4ubv((GLubyte *) & CTX::instance()->color.mesh.normals);
for(unsigned int i = 0; i < elements.size(); i++){
- MElement *ele = elements[i];
+ MElement *ele = (MElement*) elements[i];
if(!isElementVisible(ele)) continue;
SVector3 n = ele->getFace(0).normal();
for(int j = 0; j < 3; j++)
@@ -311,7 +312,7 @@ static void drawVerticesPerElement(drawContext *ctx, GEntity *e,
std::vector<T*> &elements)
{
for(unsigned int i = 0; i < elements.size(); i++){
- MElement *ele = elements[i];
+ MElement *ele = (MElement*) elements[i];
for(int j = 0; j < ele->getNumVertices(); j++){
MVertex *v = ele->getVertex(j);
if(isElementVisible(ele) && v->getVisibility()){
@@ -345,7 +346,7 @@ static void drawBarycentricDual(std::vector<T*> &elements)
gl2psEnable(GL2PS_LINE_STIPPLE);
glBegin(GL_LINES);
for(unsigned int i = 0; i < elements.size(); i++){
- MElement *ele = elements[i];
+ MElement *ele = (MElement*) elements[i];
if(!isElementVisible(ele)) continue;
SPoint3 pc = ele->barycenter();
if(ele->getDim() == 2){
@@ -433,7 +434,7 @@ 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];
+ MElement *ele = (MElement*) elements[i];
SPoint3 pc(0., 0., 0.);
if(CTX::instance()->mesh.explode != 1.) pc = ele->barycenter();
for(int j = 0; j < ele->getNumFacesRep(); j++){
@@ -682,6 +683,7 @@ class initSmoothNormalsGFace {
{
addSmoothNormals(f, f->triangles);
addSmoothNormals(f, f->quadrangles);
+ addSmoothNormals(f, f->polygons);
}
};
@@ -692,7 +694,7 @@ class initMeshGFace {
{
int num = 0;
if(CTX::instance()->mesh.surfacesEdges){
- num += (3 * f->triangles.size() + 4 * f->quadrangles.size()) / 2;
+ 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;
}
@@ -702,7 +704,7 @@ class initMeshGFace {
{
int num = 0;
if(CTX::instance()->mesh.surfacesFaces){
- num += (f->triangles.size() + 2 * f->quadrangles.size());
+ num += (f->triangles.size() + 2 * f->quadrangles.size() + 2 * f->polygons.size());
if(_curved) num *= 4;
}
return num + 100;
@@ -773,17 +775,20 @@ class drawMeshGFace {
drawVerticesPerElement(_ctx, f, f->triangles);
if(CTX::instance()->mesh.quadrangles)
drawVerticesPerElement(_ctx, f, f->quadrangles);
+ drawVerticesPerElement(_ctx, f, f->polygons);
}
}
if(CTX::instance()->mesh.normals) {
if(CTX::instance()->mesh.triangles) drawNormals(_ctx, f->triangles);
if(CTX::instance()->mesh.quadrangles) drawNormals(_ctx, f->quadrangles);
+ drawNormals(_ctx, f->polygons);
}
if(CTX::instance()->mesh.dual) {
if(CTX::instance()->mesh.triangles) drawBarycentricDual(f->triangles);
if(CTX::instance()->mesh.quadrangles) drawBarycentricDual(f->quadrangles);
+ drawBarycentricDual(f->polygons);
}
else if(CTX::instance()->mesh.voronoi) {
if(CTX::instance()->mesh.triangles) drawVoronoiDual(f->triangles);
@@ -817,8 +822,11 @@ class initMeshGRegion {
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()) / 4;
+ 18 * r->prisms.size() + 16 * r->pyramids.size() + numLP) / 4;
num = _estimateIfClipped(num);
if(CTX::instance()->mesh.explode != 1.) num *= 4;
if(_curved) num *= 2;
@@ -829,8 +837,11 @@ class initMeshGRegion {
{
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()) / 2;
+ 8 * r->prisms.size() + 6 * r->pyramids.size() + numFP) / 2;
num = _estimateIfClipped(num);
if(CTX::instance()->mesh.explode != 1.) num *= 2;
if(_curved) num *= 4;
@@ -920,6 +931,7 @@ class drawMeshGRegion {
if(CTX::instance()->mesh.hexahedra) drawVerticesPerElement(_ctx, r, r->hexahedra);
if(CTX::instance()->mesh.prisms) drawVerticesPerElement(_ctx, r, r->prisms);
if(CTX::instance()->mesh.pyramids) drawVerticesPerElement(_ctx, r, r->pyramids);
+ drawVerticesPerElement(_ctx, r, r->polyhedra);
}
}
@@ -928,6 +940,7 @@ class drawMeshGRegion {
if(CTX::instance()->mesh.hexahedra) drawBarycentricDual(r->hexahedra);
if(CTX::instance()->mesh.prisms) drawBarycentricDual(r->prisms);
if(CTX::instance()->mesh.pyramids) drawBarycentricDual(r->pyramids);
+ drawBarycentricDual(r->polyhedra);
}
if(select) {
diff --git a/contrib/DiscreteIntegration/DILevelset.cpp b/contrib/DiscreteIntegration/DILevelset.cpp
index 5060ce5060..03a4c5bd0d 100644
--- a/contrib/DiscreteIntegration/DILevelset.cpp
+++ b/contrib/DiscreteIntegration/DILevelset.cpp
@@ -5,31 +5,31 @@
#include <queue>
#include <stack>
-inline double det3 (double d11, double d12, double d13, double d21, double d22, double d23, double d31, double d32, double d33) {
- return d11*(d22*d33-d23*d32)-d21*(d12*d33-d13*d32)+d31*(d12*d23-d13*d22);
+inline double det3(double d11, double d12, double d13, double d21, double d22, double d23, double d31, double d32, double d33) {
+ return d11 * (d22 * d33 - d23 * d32) - d21 * (d12 * d33 - d13 * d32) + d31 * (d12 * d23 - d13 * d22);
}
-inline void norm (const double *vec, double *norm) {
- double n = sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
- norm[0] = vec[0]/n;
- norm[1] = vec[1]/n;
- norm[2] = vec[2]/n;
+inline void norm(const double *vec, double *norm) {
+ double n = sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
+ norm[0] = vec[0] / n;
+ norm[1] = vec[1] / n;
+ norm[2] = vec[2] / n;
}
-inline void cross (const double *pt0, const double *pt1, const double *pt2, double *cross) {
- double v1[3] = {pt1[0]-pt0[0],pt1[1]-pt0[1],pt1[2]-pt0[2]};
- double v2[3] = {pt2[0]-pt0[0],pt2[1]-pt0[1],pt2[2]-pt0[2]};
- cross[0] = v1[1]*v2[2]-v2[1]*v1[2];
- cross[1] = v2[0]*v1[2]-v1[0]*v2[2];
- cross[2] = v1[0]*v2[1]-v2[0]*v1[1];
+inline void cross(const double *pt0, const double *pt1, const double *pt2, double *cross) {
+ double v1[3] = {pt1[0] - pt0[0], pt1[1] - pt0[1], pt1[2] - pt0[2]};
+ double v2[3] = {pt2[0] - pt0[0], pt2[1] - pt0[1], pt2[2] - pt0[2]};
+ cross[0] = v1[1] * v2[2] - v2[1] * v1[2];
+ cross[1] = v2[0] * v1[2] - v1[0] * v2[2];
+ cross[2] = v1[0] * v2[1] - v2[0] * v1[1];
}
-inline bool isPlanar (const double *pt1, const double *pt2, const double *pt3, const double *pt4) {
- double c1[3]; cross(pt1,pt2,pt3,c1);
- double c2[3]; cross(pt1,pt2,pt4,c2);
- double n1[3]; norm(c1,n1);
- double n2[3]; norm(c2,n2);
- return (n1[0]==n2[0] && n1[1]==n2[1] && n1[2]==n2[2]);
+inline bool isPlanar(const double *pt1, const double *pt2, const double *pt3, const double *pt4) {
+ double c1[3]; cross(pt1, pt2, pt3, c1);
+ double c2[3]; cross(pt1, pt2, pt4, c2);
+ double n1[3]; norm(c1, n1);
+ double n2[3]; norm(c2, n2);
+ return (n1[0] == n2[0] && n1[1] == n2[1] && n1[2] == n2[2]);
}
// extrude a list of the primitive levelsets with a "Level-order traversal sequence"
@@ -44,7 +44,7 @@ void gLevelset::getPrimitives(std::vector<const gLevelset *> &gLsPrimitives) con
gLsPrimitives.push_back(p);
Q.pop();
if(!pp.empty()){
- for(int i=0; i<(int)pp.size();i++)
+ for(int i = 0; i < (int)pp.size(); i++)
Q.push(pp[i]);
}
}
@@ -53,7 +53,7 @@ void gLevelset::getPrimitives(std::vector<const gLevelset *> &gLsPrimitives) con
// return a list with the levelsets in a "Reverse Polish Notation"
void gLevelset::getRPN(std::vector<const gLevelset *> &gLsRPN) const {
std::stack<const gLevelset *> S;
- std::stack<const gLevelset *> Sc; // levelset checked
+ std::stack<const gLevelset *> Sc; // levelset checked
S.push(this);
while(!S.empty()){
const gLevelset *p = S.top();
@@ -64,11 +64,11 @@ void gLevelset::getRPN(std::vector<const gLevelset *> &gLsRPN) const {
S.pop();
}
else {
- if(Sc.empty() || p!=Sc.top()) {
- for(int i=1;i<(int)pp.size();i++) Sc.push(p);
- for(int i=(int)pp.size()-1; i>=0; i--) {
+ if(Sc.empty() || p != Sc.top()) {
+ for(int i = 1; i < (int)pp.size(); i++) Sc.push(p);
+ for(int i = (int)pp.size() - 1; i >= 0; i--) {
S.push(pp[i]);
- if(i>1) S.push(p);
+ if(i > 1) S.push(p);
}
}
else { // levelset has been checked
@@ -80,272 +80,280 @@ void gLevelset::getRPN(std::vector<const gLevelset *> &gLsRPN) const {
}
}
-gLevelsetPlane::gLevelsetPlane (const double * pt, const double *norm, int &tag) : gLevelsetPrimitive(tag) {
+gLevelsetPlane::gLevelsetPlane(const double * pt, const double *norm, int &tag) : gLevelsetPrimitive(tag) {
a = norm[0];
b = norm[1];
c = norm[2];
- d = -a*pt[0]-b*pt[1]-c*pt[2];
+ d = -a * pt[0] - b * pt[1] - c * pt[2];
}
-gLevelsetPlane::gLevelsetPlane (const double * pt1, const double *pt2, const double *pt3, int &tag) : gLevelsetPrimitive(tag) {
- a = det3(1,pt1[1],pt1[2],1,pt2[1],pt2[2],1,pt3[1],pt3[2]);
- b = det3(pt1[0],1,pt1[2],pt2[0],1,pt2[2],pt3[0],1,pt3[2]);
- c = det3(pt1[0],pt1[1],1,pt2[0],pt2[1],1,pt3[0],pt3[1],1);
- d = -det3(pt1[0],pt1[1],pt1[2],pt2[0],pt2[1],pt2[2],pt3[0],pt3[1],pt3[2]);
- //printf("a=%g,b=%g,c=%g,d=%g\n",a,b,c,d);
+gLevelsetPlane::gLevelsetPlane(const double * pt1, const double *pt2, const double *pt3, int &tag) : gLevelsetPrimitive(tag) {
+ a = det3(1., pt1[1], pt1[2], 1., pt2[1], pt2[2], 1., pt3[1], pt3[2]);
+ b = det3(pt1[0], 1., pt1[2], pt2[0], 1., pt2[2], pt3[0], 1., pt3[2]);
+ c = det3(pt1[0], pt1[1], 1., pt2[0], pt2[1], 1., pt3[0], pt3[1], 1.);
+ d = -det3(pt1[0], pt1[1], pt1[2], pt2[0], pt2[1], pt2[2], pt3[0], pt3[1], pt3[2]);
}
/*
assume a quadric
x^T A x + b^T x + c = 0
- typically , we add a rotation x -> R x
+ typically, we add a rotation x -> R x
x^T (R^T A R) x + b^T R x + c = 0
and a translation x -> x - t
x^T A x + [b^T - 2 A t] x + [c - b^T t + t^T A t ] = 0
*/
-void gLevelsetQuadric::Ax (const double x[3], double res[3], double fact){
- for (int i=0;i<3;i++){
- res[i]=0.;
- for (int j=0;j<3;j++){
- res[i] += A[i][j] * x [j] * fact;
+void gLevelsetQuadric::Ax(const double x[3], double res[3], double fact){
+ for(int i = 0; i < 3; i++){
+ res[i] = 0.;
+ for(int j = 0; j < 3; j++){
+ res[i] += A[i][j] * x[j] * fact;
}
}
}
-void gLevelsetQuadric::xAx (const double x[3], double &res, double fact){
- res= fact * (A[0][0] * x[0]*x[0] + A[1][1] * x[1]*x[1] + A[2][2] * x[2]*x[2]
- + A[1][0] * x[1]*x[0] *2 + A[2][0] * x[2]*x[0] *2 + A[1][2] * x[1]*x[2] *2);
+void gLevelsetQuadric::xAx(const double x[3], double &res, double fact){
+ res= fact * (A[0][0] * x[0] * x[0] + A[1][1] * x[1] * x[1] + A[2][2] * x[2] * x[2]
+ + A[1][0] * x[1] * x[0] * 2. + A[2][0] * x[2] * x[0] * 2. + A[1][2] * x[1] * x[2] * 2.);
}
-void gLevelsetQuadric::translate (const double transl[3]){
+void gLevelsetQuadric::translate(const double transl[3]){
double b;
- xAx (transl, b, 1.0);
- C += (-B[0]*transl[0] - B[1]*transl[1] - B[2]*transl[2] + b);
+ xAx(transl, b, 1.0);
+ C += (-B[0] * transl[0] - B[1] * transl[1] - B[2] * transl[2] + b);
double x[3];
- Ax (transl , x , 2.0);
+ Ax(transl, x, 2.0);
B[0] += -x[0];
B[1] += -x[1];
B[2] += -x[2];
}
-void gLevelsetQuadric::rotate (const double rotate[3][3]){
- double a11=0.,a12=0.,a13=0.,a22=0.,a23=0.,a33=0.,b1=0.,b2=0.,b3=0.;
- for(int i=0;i<3;i++){
- b1 += B[i]*rotate[i][0];
- b2 += B[i]*rotate[i][1];
- b3 += B[i]*rotate[i][2];
- for(int j=0;j<3;j++){
- a11 += rotate[i][0]*A[i][j]*rotate[j][0];
- a12 += rotate[i][0]*A[i][j]*rotate[j][1];
- a13 += rotate[i][0]*A[i][j]*rotate[j][2];
- a22 += rotate[i][1]*A[i][j]*rotate[j][1];
- a23 += rotate[i][1]*A[i][j]*rotate[j][2];
- a33 += rotate[i][2]*A[i][j]*rotate[j][2];
+void gLevelsetQuadric::rotate(const double rotate[3][3]){
+ double a11 = 0., a12 = 0., a13 = 0., a22 = 0., a23 = 0., a33 = 0., b1 = 0., b2 = 0., b3 = 0.;
+ for(int i = 0; i < 3; i++){
+ b1 += B[i] * rotate[i][0];
+ b2 += B[i] * rotate[i][1];
+ b3 += B[i] * rotate[i][2];
+ for(int j = 0; j < 3; j++){
+ a11 += rotate[i][0] * A[i][j] * rotate[j][0];
+ a12 += rotate[i][0] * A[i][j] * rotate[j][1];
+ a13 += rotate[i][0] * A[i][j] * rotate[j][2];
+ a22 += rotate[i][1] * A[i][j] * rotate[j][1];
+ a23 += rotate[i][1] * A[i][j] * rotate[j][2];
+ a33 += rotate[i][2] * A[i][j] * rotate[j][2];
}
}
- A[0][0]=a11;
- A[0][1]=A[1][0]=a12;
- A[0][2]=A[2][0]=a13;
- A[1][1]=a22;
- A[1][2]=A[2][1]=a23;
- A[2][2]=a33;
- B[0]=b1;
- B[1]=b2;
- B[2]=b3;
+ A[0][0] = a11;
+ A[0][1] = A[1][0] = a12;
+ A[0][2] = A[2][0] = a13;
+ A[1][1] = a22;
+ A[1][2] = A[2][1] = a23;
+ A[2][2] = a33;
+ B[0] = b1;
+ B[1] = b2;
+ B[2] = b3;
}
// computes the rotation matrix of the rotation of the vector (0,0,1) to dir
-void gLevelsetQuadric::computeRotationMatrix (const double dir[3], double t[3][3]){
- double norm = sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
- double n[3]={1.,0.,0.};
- double ct=1.,st=0.;
- if(norm !=0.) {
- double theta = atan(norm/dir[2]);
- n[0]=dir[1]/norm; n[1]=-dir[0]/norm;
+void gLevelsetQuadric::computeRotationMatrix(const double dir[3], double t[3][3]){
+ double norm = sqrt(dir[0] * dir[0] + dir[1] * dir[1]);
+ double n[3] = {1., 0., 0.};
+ double ct = 1., st = 0.;
+ if(norm != 0.) {
+ double theta = atan(norm / dir[2]);
+ n[0] = dir[1] / norm; n[1] = -dir[0] / norm;
ct = cos(theta);
st = sin(theta);
}
- t[0][0]=n[0]*n[0]*(1-ct)+ct;
- t[0][1]=n[0]*n[1]*(1-ct)-n[2]*st;
- t[0][2]=n[0]*n[2]*(1-ct)+n[1]*st;
- t[1][0]=n[1]*n[0]*(1-ct)+n[2]*st;
- t[1][1]=n[1]*n[1]*(1-ct)+ct;
- t[1][2]=n[1]*n[2]*(1-ct)-n[0]*st;
- t[2][0]=n[2]*n[0]*(1-ct)-n[1]*st;
- t[2][1]=n[2]*n[1]*(1-ct)+n[0]*st;
- t[2][2]=n[2]*n[2]*(1-ct)+ct;
+ t[0][0] = n[0] * n[0] * (1. - ct) + ct;
+ t[0][1] = n[0] * n[1] * (1. - ct) - n[2] * st;
+ t[0][2] = n[0] * n[2] * (1. - ct) + n[1] * st;
+ t[1][0] = n[1] * n[0] * (1. - ct) + n[2] * st;
+ t[1][1] = n[1] * n[1] * (1. - ct) + ct;
+ t[1][2] = n[1] * n[2] * (1. - ct) - n[0] * st;
+ t[2][0] = n[2] * n[0] * (1. - ct) - n[1] * st;
+ t[2][1] = n[2] * n[1] * (1. - ct) + n[0] * st;
+ t[2][2] = n[2] * n[2] * (1. - ct) + ct;
}
-void gLevelsetQuadric::computeRotationMatrix (const double dir1[3], const double dir2[3] , double t[3][3]){
- double norm = sqrt((dir1[1]*dir2[2]-dir1[2]*dir2[1])*(dir1[1]*dir2[2]-dir1[2]*dir2[1])
- + (dir1[2]*dir2[0]-dir1[0]*dir2[2])*(dir1[2]*dir2[0]-dir1[0]*dir2[2])
- + (dir1[0]*dir2[1]-dir1[1]*dir2[0])*(dir1[0]*dir2[1]-dir1[1]*dir2[0]));
- double n[3]={1.,0.,0.};
- double ct=1.,st=0.;
- if(norm!=0.) {
- st = norm/((dir1[0]*dir1[0]+dir1[1]*dir1[1]+dir1[2]*dir1[2])*(dir2[0]*dir2[0]+dir2[1]*dir2[1]+dir2[2]*dir2[2]));
- n[0]=(dir1[1]*dir2[2]-dir1[2]*dir2[1])/norm;
- n[1]=(dir1[2]*dir2[0]-dir1[0]*dir2[2])/norm;
- n[2]=(dir1[0]*dir2[1]-dir1[1]*dir2[0])/norm;
+void gLevelsetQuadric::computeRotationMatrix(const double dir1[3], const double dir2[3], double t[3][3]){
+ double norm = sqrt((dir1[1] * dir2[2] - dir1[2] * dir2[1]) * (dir1[1] * dir2[2] - dir1[2] * dir2[1])
+ + (dir1[2] * dir2[0] - dir1[0] * dir2[2]) * (dir1[2] * dir2[0] - dir1[0] * dir2[2])
+ + (dir1[0] * dir2[1] - dir1[1] * dir2[0]) * (dir1[0] * dir2[1] - dir1[1] * dir2[0]));
+ double n[3] = {1., 0., 0.};
+ double ct = 1., st = 0.;
+ if(norm != 0.) {
+ st = norm / ((dir1[0] * dir1[0] + dir1[1] * dir1[1] + dir1[2] * dir1[2]) *
+ (dir2[0] * dir2[0] + dir2[1] * dir2[1] + dir2[2] * dir2[2]));
+ n[0] = (dir1[1] * dir2[2] - dir1[2] * dir2[1]) / norm;
+ n[1] = (dir1[2] * dir2[0] - dir1[0] * dir2[2]) / norm;
+ n[2] = (dir1[0] * dir2[1] - dir1[1] * dir2[0]) / norm;
ct = cos(asin(st));
}
- t[0][0]=n[0]*n[0]*(1-ct)+ct;
- t[0][1]=n[0]*n[1]*(1-ct)-n[2]*st;
- t[0][2]=n[0]*n[2]*(1-ct)+n[1]*st;
- t[1][0]=n[1]*n[0]*(1-ct)+n[2]*st;
- t[1][1]=n[1]*n[1]*(1-ct)+ct;
- t[1][2]=n[1]*n[2]*(1-ct)-n[0]*st;
- t[2][0]=n[2]*n[0]*(1-ct)-n[1]*st;
- t[2][1]=n[2]*n[1]*(1-ct)+n[0]*st;
- t[2][2]=n[2]*n[2]*(1-ct)+ct;
+ t[0][0] = n[0] * n[0] * (1. - ct) + ct;
+ t[0][1] = n[0] * n[1] * (1. - ct) - n[2] * st;
+ t[0][2] = n[0] * n[2] * (1. - ct) + n[1] * st;
+ t[1][0] = n[1] * n[0] * (1. - ct) + n[2] * st;
+ t[1][1] = n[1] * n[1] * (1. - ct) + ct;
+ t[1][2] = n[1] * n[2] * (1. - ct) - n[0] * st;
+ t[2][0] = n[2] * n[0] * (1. - ct) - n[1] * st;
+ t[2][1] = n[2] * n[1] * (1. - ct) + n[0] * st;
+ t[2][2] = n[2] * n[2] * (1. - ct) + ct;
}
void gLevelsetQuadric::init(){
- for (int i=0;i<3;i++){
+ for(int i = 0; i < 3; i++){
B[i] = 0.;
- for (int j=0;j<3;j++)A[i][j]=0.;
+ for(int j = 0; j < 3; j++) A[i][j] = 0.;
}
C = 0.;
}
-double gLevelsetQuadric::operator () (const double &x, const double &y, const double &z) const{
- return (A[0][0]*x*x + 2*A[0][1]*x*y + 2*A[0][2]*x*z + A[1][1]*y*y + 2*A[1][2]*y*z + A[2][2]*z*z + B[0]*x + B[1]*y + B[2]*z + C);
+double gLevelsetQuadric::operator()(const double &x, const double &y, const double &z) const{
+ return(A[0][0] * x * x + 2. * A[0][1] * x * y + 2. * A[0][2] * x * z + A[1][1] * y * y
+ + 2. * A[1][2] * y * z + A[2][2] * z * z + B[0] * x + B[1] * y + B[2] * z + C);
}
-gLevelsetGenCylinder :: gLevelsetGenCylinder (const double * pt ,const double *dir ,const double &R, int &tag) : gLevelsetQuadric(tag) {
+gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir, const double &R,
+ int &tag) : gLevelsetQuadric(tag) {
A[0][0] = 1.;
A[1][1] = 1.;
- C = - R*R;
+ C = - R * R;
double rot[3][3];
- computeRotationMatrix(dir,rot);
+ computeRotationMatrix(dir, rot);
rotate(rot);
- translate (pt);
+ translate(pt);
}
-gLevelsetEllipsoid::gLevelsetEllipsoid (const double *pt, const double *dir, const double &a, const double &b, const double &c, int &tag) : gLevelsetQuadric(tag) {
- A[0][0] = 1./(a*a);
- A[1][1] = 1./(b*b);
- A[2][2] = 1./(c*c);
+gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir, const double &a,
+ const double &b, const double &c, int &tag) : gLevelsetQuadric(tag) {
+ A[0][0] = 1. / (a * a);
+ A[1][1] = 1. / (b * b);
+ A[2][2] = 1. / (c * c);
C = -1.;
- double rot [3][3];
- computeRotationMatrix(dir,rot);
+ double rot[3][3];
+ computeRotationMatrix(dir, rot);
rotate(rot);
translate(pt);
}
-gLevelsetCone::gLevelsetCone (const double * pt, const double * dir, const double &angle, int &tag) : gLevelsetQuadric(tag) {
+gLevelsetCone::gLevelsetCone(const double *pt, const double *dir, const double &angle, int &tag) : gLevelsetQuadric(tag) {
A[0][0] = 1.;
A[1][1] = 1.;
- A[2][2] = -tan(angle)*tan(angle);
- double rot [3][3];
- computeRotationMatrix(dir,rot);
+ A[2][2] = -tan(angle) * tan(angle);
+ double rot[3][3];
+ computeRotationMatrix(dir, rot);
rotate(rot);
translate(pt);
}
-gLevelsetGeneralQuadric::gLevelsetGeneralQuadric (const double *pt, const double *dir, const double &x2, const double &y2, const double &z2,
- const double &z, const double &c, int &tag) : gLevelsetQuadric(tag) {
+gLevelsetGeneralQuadric::gLevelsetGeneralQuadric(const double *pt, const double *dir, const double &x2, const double &y2, const double &z2,
+ const double &z, const double &c, int &tag) : gLevelsetQuadric(tag) {
A[0][0] = x2;
A[1][1] = y2;
A[2][2] = z2;
B[2] = z;
- C=c;
- double rot [3][3];
- computeRotationMatrix(dir,rot);
+ C = c;
+ double rot[3][3];
+ computeRotationMatrix(dir, rot);
rotate(rot);
translate(pt);
}
gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *dir2, const double *dir3,
const double &a, const double &b, const double &c, int &tag) {
- double dir1m[3]={-dir1[0],-dir1[1],-dir1[2]};
- double dir2m[3]={-dir2[0],-dir2[1],-dir2[2]};
- double dir3m[3]={-dir3[0],-dir3[1],-dir3[2]};
- double n1[3]; norm(dir1,n1);
- double n2[3]; norm(dir2,n2);
- double n3[3]; norm(dir3,n3);
- double pt2[3]={pt[0]+a*n1[0]+b*n2[0]+c*n3[0],pt[1]+a*n1[1]+b*n2[1]+c*n3[1],pt[2]+a*n1[2]+b*n2[2]+c*n3[2]};
+ double dir1m[3] = {-dir1[0], -dir1[1], -dir1[2]};
+ double dir2m[3] = {-dir2[0], -dir2[1], -dir2[2]};
+ double dir3m[3] = {-dir3[0], -dir3[1], -dir3[2]};
+ double n1[3]; norm(dir1, n1);
+ double n2[3]; norm(dir2, n2);
+ double n3[3]; norm(dir3, n3);
+ double pt2[3] = {pt[0] + a * n1[0] + b * n2[0] + c * n3[0], pt[1] + a * n1[1] + b * n2[1] + c * n3[1],
+ pt[2] + a * n1[2] + b * n2[2] + c * n3[2]};
std::vector<const gLevelset *> p;
- p.push_back(new gLevelsetPlane (pt2,dir3,tag));
- p.push_back(new gLevelsetPlane (pt,dir3m,tag));
- p.push_back(new gLevelsetPlane (pt,dir2m,tag));
- p.push_back(new gLevelsetPlane (pt2,dir2,tag));
- p.push_back(new gLevelsetPlane (pt2,dir1,tag));
- p.push_back(new gLevelsetPlane (pt,dir1m,tag));
+ p.push_back(new gLevelsetPlane(pt2, dir3, tag));
+ p.push_back(new gLevelsetPlane(pt, dir3m, tag));
+ p.push_back(new gLevelsetPlane(pt, dir2m, tag));
+ p.push_back(new gLevelsetPlane(pt2, dir2, tag));
+ p.push_back(new gLevelsetPlane(pt2, dir1, tag));
+ p.push_back(new gLevelsetPlane(pt, dir1m, tag));
Ls = new gLevelsetIntersection(p);
}
gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *pt3, const double *pt4,
const double *pt5, const double *pt6, const double *pt7, const double *pt8, int &tag) {
- if(!isPlanar(pt1,pt2,pt3,pt4) || !isPlanar(pt5,pt6,pt7,pt8) || !isPlanar(pt1,pt2,pt5,pt6) ||
- !isPlanar(pt3,pt4,pt7,pt8) || !isPlanar(pt1,pt4,pt5,pt8) || !isPlanar(pt2,pt3,pt6,pt7))
- printf("WARNING : faces of the box are not planar! %d,%d,%d,%d,%d,%d\n",
- isPlanar(pt1,pt2,pt3,pt4), isPlanar(pt5,pt6,pt7,pt8), isPlanar(pt1,pt2,pt5,pt6),
- isPlanar(pt3,pt4,pt7,pt8), isPlanar(pt1,pt4,pt5,pt8), isPlanar(pt2,pt3,pt6,pt7));
+ if(!isPlanar(pt1, pt2, pt3, pt4) || !isPlanar(pt5, pt6, pt7, pt8) || !isPlanar(pt1, pt2, pt5, pt6) ||
+ !isPlanar(pt3, pt4, pt7, pt8) || !isPlanar(pt1, pt4, pt5, pt8) || !isPlanar(pt2, pt3, pt6, pt7))
+ printf("WARNING : faces of the box are not planar! %d, %d, %d, %d, %d, %d\n",
+ isPlanar(pt1, pt2, pt3, pt4), isPlanar(pt5, pt6, pt7, pt8), isPlanar(pt1, pt2, pt5, pt6),
+ isPlanar(pt3, pt4, pt7, pt8), isPlanar(pt1, pt4, pt5, pt8), isPlanar(pt2, pt3, pt6, pt7));
std::vector<const gLevelset *> p;
- p.push_back(new gLevelsetPlane (pt5,pt6,pt8,tag));
- p.push_back(new gLevelsetPlane (pt1,pt4,pt2,tag));
- p.push_back(new gLevelsetPlane (pt1,pt2,pt5,tag));
- p.push_back(new gLevelsetPlane (pt3,pt4,pt7,tag));
- p.push_back(new gLevelsetPlane (pt2,pt3,pt6,tag));
- p.push_back(new gLevelsetPlane (pt1,pt5,pt4,tag));
+ p.push_back(new gLevelsetPlane(pt5, pt6, pt8, tag));
+ p.push_back(new gLevelsetPlane(pt1, pt4, pt2, tag));
+ p.push_back(new gLevelsetPlane(pt1, pt2, pt5, tag));
+ p.push_back(new gLevelsetPlane(pt3, pt4, pt7, tag));
+ p.push_back(new gLevelsetPlane(pt2, pt3, pt6, tag));
+ p.push_back(new gLevelsetPlane(pt1, pt5, pt4, tag));
Ls = new gLevelsetIntersection(p);
}
-gLevelsetCylinder :: gLevelsetCylinder (const double * pt ,const double *dir ,const double &R, const double &H, int &tag) {
- double dir2[3] = {-dir[0],-dir[1],-dir[2]};
- double n[3]; norm(dir,n);
- double pt2[3] = {pt[0]+H*n[0],pt[1]+H*n[1],pt[2]+H*n[2]};
+gLevelsetCylinder::gLevelsetCylinder(const double *pt, const double *dir, const double &R, const double &H, int &tag) {
+ double dir2[3] = {-dir[0], -dir[1], -dir[2]};
+ double n[3]; norm(dir, n);
+ double pt2[3] = {pt[0] + H * n[0], pt[1] + H * n[1], pt[2] + H * n[2]};
std::vector<const gLevelset *> p;
- p.push_back(new gLevelsetGenCylinder (pt,dir,R,tag));
- p.push_back(new gLevelsetPlane (pt,dir2,tag));
- p.push_back(new gLevelsetPlane (pt2,dir,tag));
+ p.push_back(new gLevelsetGenCylinder(pt, dir, R, tag));
+ p.push_back(new gLevelsetPlane(pt, dir2, tag));
+ p.push_back(new gLevelsetPlane(pt2, dir, tag));
Ls = new gLevelsetIntersection(p);
}
-gLevelsetCylinder :: gLevelsetCylinder (const double * pt ,const double *dir ,const double &R, const double &r, const double &H, int &tag) {
- double dir2[3] = {-dir[0],-dir[1],-dir[2]};
- double n[3]; norm(dir,n);
- double pt2[3] = {pt[0]+H*n[0],pt[1]+H*n[1],pt[2]+H*n[2]};
+gLevelsetCylinder::gLevelsetCylinder(const double * pt, const double *dir, const double &R, const double &r, const double &H, int &tag) {
+ double dir2[3] = {-dir[0], -dir[1], -dir[2]};
+ double n[3]; norm(dir, n);
+ double pt2[3] = {pt[0] + H * n[0], pt[1] + H * n[1], pt[2] + H * n[2]};
std::vector<const gLevelset *> p1;
- p1.push_back(new gLevelsetGenCylinder (pt,dir,R,tag));
- p1.push_back(new gLevelsetPlane (pt,dir2,tag));
- p1.push_back(new gLevelsetPlane (pt2,dir,tag));
+ p1.push_back(new gLevelsetGenCylinder(pt, dir, R, tag));
+ p1.push_back(new gLevelsetPlane(pt, dir2, tag));
+ p1.push_back(new gLevelsetPlane(pt2, dir, tag));
std::vector<const gLevelset *> p2;
- p2.push_back(new gLevelsetIntersection (p1));
- p2.push_back(new gLevelsetGenCylinder (pt,dir,r,tag));
+ p2.push_back(new gLevelsetIntersection(p1));
+ p2.push_back(new gLevelsetGenCylinder(pt, dir, r, tag));
Ls = new gLevelsetCut(p2);
}
-gLevelsetConrod::gLevelsetConrod (const double *pt, const double *dir1, const double *dir2,
- const double &H1, const double &H2, const double &H3,
- const double &R1, const double &r1, const double &R2, const double &r2,
- const double &L1, const double &L2, const double &E, int &tag) {
- double n1[3]; norm(dir1,n1);
- double n2[3]; norm(dir2,n2);
- double pt1[3] = {pt[0]-n2[0]*H1/2,pt[1]-n2[1]*H1/2,pt[2]-n2[2]*H1/2};
- double pt2[3] = {pt[0]+n1[0]*E-n2[0]*H2/2,pt[1]+n1[1]*E-n2[1]*H2/2,pt[2]+n1[2]*E-n2[2]*H2/2};
- double dir3[3]; cross(pt1,pt2,pt,dir3);
- double n3[3]; norm(dir3,n3);
- double pt31[3] = {pt[0]-n2[0]*H3/2+n3[0]*L1/2,pt[1]-n2[1]*H3/2+n3[1]*L1/2,pt[2]-n2[2]*H3/2+n3[2]*L1/2};
- double pt32[3] = {pt31[0]-n3[0]*L1,pt31[1]-n3[1]*L1,pt31[2]-n3[2]*L1};
- double pt33[3] = {pt32[0]+n2[0]*H3,pt32[1]+n2[1]*H3,pt32[2]+n2[2]*H3};
- double pt34[3] = {pt33[0]+n3[0]*L1,pt33[1]+n3[1]*L1,pt33[2]+n3[2]*L1};
- double pt35[3] = {pt[0]+n1[0]*E-n2[0]*H3/2+n3[0]*L2/2,pt[1]+n1[1]*E-n2[1]*H3/2+n3[1]*L2/2,pt[2]+n1[2]*E-n2[2]*H3/2+n3[2]*L2/2};
- double pt36[3] = {pt35[0]-n3[0]*L2,pt35[1]-n3[1]*L2,pt35[2]-n3[2]*L2};
- double pt37[3] = {pt36[0]+n2[0]*H3,pt36[1]+n2[1]*H3,pt36[2]+n2[2]*H3};
- double pt38[3] = {pt37[0]+n3[0]*L2,pt37[1]+n3[1]*L2,pt37[2]+n3[2]*L2};
+gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1, const double *dir2,
+ const double &H1, const double &H2, const double &H3,
+ const double &R1, const double &r1, const double &R2, const double &r2,
+ const double &L1, const double &L2, const double &E, int &tag) {
+ double n1[3]; norm(dir1, n1);
+ double n2[3]; norm(dir2, n2);
+ double pt1[3] = {pt[0] - n2[0] * H1 / 2., pt[1] - n2[1] * H1 / 2., pt[2] - n2[2] * H1 / 2.};
+ double pt2[3] = {pt[0] + n1[0] * E - n2[0] * H2 / 2., pt[1] + n1[1] * E - n2[1] * H2 / 2.,
+ pt[2] + n1[2] * E - n2[2] * H2 / 2.};
+ double dir3[3]; cross(pt1, pt2, pt, dir3);
+ double n3[3]; norm(dir3, n3);
+ double pt31[3] = {pt[0] - n2[0] * H3 / 2. + n3[0] * L1 / 2., pt[1] - n2[1] * H3 / 2. + n3[1] * L1 / 2.,
+ pt[2] - n2[2] * H3 / 2. + n3[2] * L1 / 2.};
+ double pt32[3] = {pt31[0] - n3[0] * L1, pt31[1] - n3[1] * L1, pt31[2] - n3[2] * L1};
+ double pt33[3] = {pt32[0] + n2[0] * H3, pt32[1] + n2[1] * H3, pt32[2] + n2[2] * H3};
+ double pt34[3] = {pt33[0] + n3[0] * L1, pt33[1] + n3[1] * L1, pt33[2] + n3[2] * L1};
+ double pt35[3] = {pt[0] + n1[0] * E - n2[0] * H3 / 2. + n3[0] * L2 / 2.,
+ pt[1] + n1[1] * E - n2[1] * H3 / 2. + n3[1] * L2 / 2.,
+ pt[2] + n1[2] * E - n2[2] * H3 / 2. + n3[2] * L2 / 2.};
+ double pt36[3] = {pt35[0] - n3[0] * L2, pt35[1] - n3[1] * L2, pt35[2] - n3[2] * L2};
+ double pt37[3] = {pt36[0] + n2[0] * H3, pt36[1] + n2[1] * H3, pt36[2] + n2[2] * H3};
+ double pt38[3] = {pt37[0] + n3[0] * L2, pt37[1] + n3[1] * L2, pt37[2] + n3[2] * L2};
std::vector<const gLevelset *> p1;
- p1.push_back(new gLevelsetBox (pt31,pt32,pt33,pt34,pt35,pt36,pt37,pt38,tag));
- p1.push_back(new gLevelsetCylinder (pt1,dir2,R1,H1,tag));
- p1.push_back(new gLevelsetCylinder (pt2,dir2,R2,H2,tag));
+ p1.push_back(new gLevelsetBox(pt31, pt32, pt33, pt34, pt35, pt36, pt37, pt38, tag));
+ p1.push_back(new gLevelsetCylinder(pt1, dir2, R1, H1, tag));
+ p1.push_back(new gLevelsetCylinder(pt2, dir2, R2, H2, tag));
std::vector<const gLevelset *> p2;
- p2.push_back(new gLevelsetUnion (p1));
- p2.push_back(new gLevelsetGenCylinder (pt1,dir2,r1,tag));
- p2.push_back(new gLevelsetGenCylinder (pt2,dir2,r2,tag));
+ p2.push_back(new gLevelsetUnion(p1));
+ p2.push_back(new gLevelsetGenCylinder(pt1, dir2, r1, tag));
+ p2.push_back(new gLevelsetGenCylinder(pt2, dir2, r2, tag));
Ls = new gLevelsetCut(p2);
}
diff --git a/contrib/DiscreteIntegration/DILevelset.h b/contrib/DiscreteIntegration/DILevelset.h
index 4bf55ea5dc..05262dc27d 100644
--- a/contrib/DiscreteIntegration/DILevelset.h
+++ b/contrib/DiscreteIntegration/DILevelset.h
@@ -32,24 +32,24 @@ protected:
public:
gLevelset() : tag_(-1) {}
virtual ~gLevelset(){}
- virtual double operator () (const double &x, const double &y, const double &z) const = 0;
+ virtual double operator() (const double &x, const double &y, const double &z) const = 0;
// inline double operator () (const SPoint3 &p) const {return this->operator()(p.x(),p.y(),p.z());}
- bool isInsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z)*insideDomain > 0.;}
- bool isOutsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z)*insideDomain < 0.;}
+ bool isInsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) * insideDomain > 0.;}
+ bool isOutsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) * insideDomain < 0.;}
bool isOnBorder (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) == 0.;}
virtual std::vector<const gLevelset *> getChildren() const = 0;
- virtual double choose (double d1, double d2) const=0;
+ virtual double choose (double d1, double d2) const = 0;
virtual int type() const = 0;
- virtual bool isPrimitive() const =0;
+ virtual bool isPrimitive() const = 0;
void setTag(int t) {
assert(isPrimitive());
- tag_=t;
+ tag_ = t;
}
virtual int getTag() const { return tag_; }
void getPrimitives(std::vector<const gLevelset *> &primitives) const;
void getRPN(std::vector<const gLevelset *> &gLsRPN) const;
double H (const double &x, const double &y, const double &z) const {
- if ( isInsideDomain(x,y,z) || isOnBorder(x,y,z)) return 1.0;
+ if (isInsideDomain(x,y,z) || isOnBorder(x,y,z)) return 1.0;
return 0.0;
}
void print() const {
@@ -69,7 +69,7 @@ public:
case INTER : printf("INTER"); break;
case LSMESH: printf("LSMESH"); break;
}
- printf(" Tag=%d\n",getTag());
+ printf(" Tag=%d\n", getTag());
}
};
@@ -80,8 +80,8 @@ class gLevelsetPrimitive : public gLevelset
public:
gLevelsetPrimitive() : gLevelset() {}
gLevelsetPrimitive(int &tag) {
- if (tag<1) {
- printf("Tag of the levelset (%d) must be greater than 0.\n",tag);
+ if (tag < 1) {
+ printf("Tag of the levelset (%d) must be greater than 0.\n", tag);
throw;
}
tag_ = tag++;
@@ -98,36 +98,36 @@ public:
class gLevelsetSphere : public gLevelsetPrimitive
{
- double xc,yc,zc,r;
+ double xc, yc, zc, r;
public:
- gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int &tag) : gLevelsetPrimitive(tag),xc(x),yc(y),zc(z),r(R) {}
+ gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int &tag) : gLevelsetPrimitive(tag), xc(x), yc(y), zc(z), r(R) {}
virtual double operator () (const double &x, const double &y, const double &z) const
- {return sqrt((xc-x)*(xc-x)+(yc-y)*(yc-y)+(zc-z)*(zc-z)) - r;}
+ {return sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) + (zc - z) * (zc - z)) - r;}
int type() const {return SPHERE;}
};
class gLevelsetPlane : public gLevelsetPrimitive
{
- double a,b,c,d;
+ double a, b, c, d;
public:
// define the plane _a*x+_b*y+_c*z+_d, with outward normal (a,b,c)
- gLevelsetPlane (const double &_a, const double &_b, const double &_c, const double &_d, int &tag) : gLevelsetPrimitive(tag),a(_a),b(_b),c(_c),d(_d) {}
+ gLevelsetPlane (const double &_a, const double &_b, const double &_c, const double &_d, int &tag) : gLevelsetPrimitive(tag), a(_a), b(_b), c(_c), d(_d) {}
// define the plane passing through the point pt and with outward normal norm
- gLevelsetPlane (const double * pt, const double *norm, int &tag);
+ gLevelsetPlane (const double *pt, const double *norm, int &tag);
// define the plane passing through the 3 points pt1,pt2,pt3 and with outward normal (pt1,pt2)x(pt1,pt3)
gLevelsetPlane (const double *pt1, const double *pt2, const double *pt3, int &tag);
// return negative value inward and positive value outward
- virtual double operator () (const double &x, const double &y, const double &z) const
- {return a*x + b*y + c*z + d;}
+ virtual double operator() (const double &x, const double &y, const double &z) const
+ {return a * x + b * y + c * z + d;}
int type() const {return PLANE;}
};
class gLevelsetQuadric : public gLevelsetPrimitive
{
protected:
- double A[3][3],B[3],C;
+ double A[3][3], B[3], C;
void translate (const double transl[3]);
- void rotate (const double rotate[3][3]);
+ void rotate (const double rotate[3][3]);
void computeRotationMatrix (const double dir[3], double t[3][3]);
void computeRotationMatrix (const double dir1[3], const double dir2[3] , double t[3][3]);
void Ax (const double x[3], double res[3], double fact=1.0);
@@ -139,27 +139,27 @@ public:
gLevelsetQuadric(int &tag) : gLevelsetPrimitive(tag) {init(); }
~gLevelsetQuadric() {}
double operator () (const double &x, const double &y, const double &z) const;
- virtual int type() const =0;
+ virtual int type() const = 0;
};
class gLevelsetGenCylinder : public gLevelsetQuadric
{
public:
- gLevelsetGenCylinder (const double * pt, const double *dir, const double &R, int &tag);
+ gLevelsetGenCylinder (const double *pt, const double *dir, const double &R, int &tag);
int type() const {return GENCYLINDER;}
};
class gLevelsetEllipsoid : public gLevelsetQuadric
{
public:
- gLevelsetEllipsoid (const double * pt, const double *dir, const double &a, const double &b, const double &c, int &tag);
+ gLevelsetEllipsoid (const double *pt, const double *dir, const double &a, const double &b, const double &c, int &tag);
int type() const {return ELLIPS;}
};
class gLevelsetCone : public gLevelsetQuadric
{
public:
- gLevelsetCone (const double * pt, const double * dir, const double &angle, int &tag);
+ gLevelsetCone (const double *pt, const double *dir, const double &angle, int &tag);
int type() const {return CONE;}
};
@@ -179,21 +179,21 @@ protected:
std::vector<const gLevelset *> children;
public:
gLevelsetTools () {}
- gLevelsetTools (std::vector<const gLevelset *> &p) {children=p;}
+ gLevelsetTools (std::vector<const gLevelset *> &p) {children = p;}
~gLevelsetTools () {
- for(int i=0;i<(int)children.size();i++)
+ for(int i = 0; i < (int)children.size(); i++)
delete children[i];
}
double operator () (const double &x, const double &y, const double &z) const {
- double d = (*children[0])(x,y,z);
- for (int i=1;i<(int)children.size();i++){
- double dt = (*children[i])(x,y,z);
- d = choose(d,dt);
+ double d = (*children[0])(x, y, z);
+ for (int i = 1; i < (int)children.size(); i++){
+ double dt = (*children[i])(x, y, z);
+ d = choose(d, dt);
}
return d;
}
std::vector<const gLevelset *> getChildren() const {return children;}
- virtual double choose (double d1, double d2) const=0;
+ virtual double choose (double d1, double d2) const = 0;
virtual int type() const = 0;
bool isPrimitive() const {return false;}
};
@@ -202,10 +202,9 @@ class gLevelsetReverse : public gLevelset
{
const gLevelset *ls;
public:
- gLevelsetReverse ( const gLevelset * p) : ls(p){
- }
+ gLevelsetReverse (const gLevelset *p) : ls(p){}
double operator () (const double &x, const double &y, const double &z) const {
- return -(*ls)(x,y,z);
+ return -(*ls)(x, y, z);
}
std::vector<const gLevelset *> getChildren() const {return ls->getChildren();}
bool isPrimitive() const {return ls->isPrimitive();}
@@ -219,9 +218,9 @@ public:
class gLevelsetCut : public gLevelsetTools
{
public:
- gLevelsetCut ( std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
+ gLevelsetCut (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
double choose (double d1, double d2) const {
- return (d1>-d2)?d1:-d2; // greater of d1 and -d2
+ return (d1 > -d2) ? d1 : -d2; // greater of d1 and -d2
}
int type() const {return CUT;}
};
@@ -230,9 +229,9 @@ public:
class gLevelsetUnion : public gLevelsetTools
{
public:
- gLevelsetUnion ( std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
+ gLevelsetUnion (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
double choose (double d1, double d2) const {
- return (d1<d2)?d1:d2; // lesser of d1 and d2
+ return (d1 < d2) ? d1 : d2; // lesser of d1 and d2
}
int type() const {return UNION;}
};
@@ -241,9 +240,9 @@ public:
class gLevelsetIntersection : public gLevelsetTools
{
public:
- gLevelsetIntersection ( std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
+ gLevelsetIntersection (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
double choose (double d1, double d2) const {
- return (d1>d2)?d1:d2; // greater of d1 and d2
+ return (d1 > d2) ? d1 : d2; // greater of d1 and d2
}
int type() const {return INTER;}
};
@@ -252,14 +251,14 @@ public:
class gLevelsetCrack : public gLevelsetTools
{
public:
- gLevelsetCrack ( std::vector<const gLevelset *> &p) {
- assert (p.size()==2);
+ gLevelsetCrack (std::vector<const gLevelset *> &p) {
+ assert (p.size() == 2);
children.push_back(p[0]);
children.push_back(new gLevelsetReverse(p[0]));
children.push_back(p[1]);
}
double choose (double d1, double d2) const {
- return (d1>d2)?d1:d2; // greater of d1 and d2
+ return (d1 > d2) ? d1 : d2; // greater of d1 and d2
}
int type() const {return CRACK;}
};
@@ -271,10 +270,10 @@ class gLevelsetImproved : public gLevelset
protected:
gLevelset *Ls;
public:
- double operator () (const double &x, const double &y, const double &z) const {return (*Ls)(x,y,z);}
+ double operator() (const double &x, const double &y, const double &z) const {return (*Ls)(x, y, z);}
std::vector<const gLevelset *> getChildren() const { return Ls->getChildren(); }
- double choose (double d1, double d2) const { return Ls->choose(d1,d2); }
- virtual int type() const =0;
+ double choose (double d1, double d2) const { return Ls->choose(d1, d2); }
+ virtual int type() const = 0;
bool isPrimitive() const {return Ls->isPrimitive();}
};
@@ -319,7 +318,7 @@ public:
// tags of the faces are : exterior face : tag+0
// plane face including pt : tag+1
// plane face opposite to pt : tag+2
- gLevelsetCylinder (const double * pt, const double *dir, const double &R, const double &H, int &tag);
+ gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &H, int &tag);
// create a cylinder : pt is the point in the middle of the cylinder base,
// dir is the direction of the cylinder axis,
// R is the outer radius of the cylinder,
@@ -329,7 +328,7 @@ public:
// plane face including pt : tag+1
// plane face opposite to pt : tag+2
// interior face : tag+3
- gLevelsetCylinder (const double * pt, const double *dir, const double &R, const double &r, const double &H, int &tag);
+ gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &r, const double &H, int &tag);
int type() const {return CYLINDER;}
};
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 1dcef683c4..4c3197eb07 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -825,7 +825,7 @@ void DI_Element::evalC (const double u, const double v, const double w, double *
int nbV = nbVert() + nbMid();
std::vector<double> s(nbV);
ev[0] = 0; ev[1] = 0; ev[2] = 0;
- getShapeFunctions (u, v, w, &s[0], order); //printf("o=%d nbV=%d s=%g,%g,%g,%g,%g,%g\n",order,nbV,s[0],s[1],s[2],s[3],s[4],s[5]);
+ getShapeFunctions (u, v, w, &s[0], order);
for(int i = 0; i < nbV; i++){
ev[0] += x(i) * s[i];
ev[1] += y(i) * s[i];
@@ -1014,15 +1014,17 @@ void DI_Line::computeIntegral() {
}
void DI_Line::getShapeFunctions (double u, double v, double w, double s[], int ord) const {
int order = (ord == -1) ? getPolynomialOrder() : ord;
+ double valm = (fabs(1. - u) < 1.e-16) ? 0. : (1. - u);
+ double valp = (fabs(1. + u) < 1.e-16) ? 0. : (1. + u);
switch (order) {
case 1 :
- s[0] = (1. - u) / 2.;
- s[1] = (1. + u) / 2.;
+ s[0] = valm / 2.;
+ s[1] = valp / 2.;
break;
case 2 :
- s[0] = u * (u - 1.) / 2.;
- s[1] = u * (u + 1.) / 2.;
- s[2] = (1. - u) * (1. + u);
+ s[0] = -u * valm / 2.;
+ s[1] = u * valp / 2.;
+ s[2] = valm * valp;
break;
default : printf("Order %d line function space not implemented ", order); print();
}
@@ -1071,19 +1073,20 @@ void DI_Triangle::computeIntegral() {
}
void DI_Triangle::getShapeFunctions (double u, double v, double w, double s[], int ord) const {
int order = (ord == -1) ? getPolynomialOrder() : ord;
+ double val1 = (fabs(1. - u - v) < 1.e-16) ? 0. : (1. - u - v);
switch (order) {
case 1 :
- s[0] = 1. - u - v;
+ s[0] = val1;
s[1] = u;
s[2] = v;
break;
case 2 :
- s[0] = (1. - u - v) * (1. - 2. * u - 2. * v);
+ s[0] = val1 * (1. - 2. * u - 2. * v);
s[1] = u * (2. * u - 1.);
s[2] = v * (2. * v - 1.);
- s[3] = 4. * u * (1. - u - v);
+ s[3] = 4. * u * val1;
s[4] = 4. * u * v;
- s[5] = 4. * v * (1. - u - v);
+ s[5] = 4. * v * val1;
break;
default :
printf("Order %d triangle function space not implemented ", order); print();
@@ -2273,7 +2276,7 @@ void DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip
assert(tt_subTriangles[t].sizeLs() == 1);
tt_subTriangles[t].computeLsTagDom(&tt, RPN);
DI_Triangle tt_subTr = tt_subTriangles[t];
- mappingEl(&tt_subTr); //tt_subTr.printls();
+ mappingEl(&tt_subTr);
tt_subTr.integrationPoints(polynomialOrderTr, &tt_subTriangles[t], &tt, RPN, ip);
subTriangles.push_back(tt_subTr);
}
--
GitLab