Skip to content
Snippets Groups Projects
Commit ad2b36e6 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

fix display problems (gbricteux)

parent 18fb762d
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/MElementCut.cpp
+ 49
− 47
View file @ ad2b36e6
......@@ -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");
......
Graphics/Makefile
+ 1
− 0
View file @ ad2b36e6
......@@ -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
......
Graphics/drawMesh.cpp
+ 21
− 8
View file @ ad2b36e6
......@@ -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) {
......
contrib/DiscreteIntegration/DILevelset.cpp
+ 199
− 191
View file @ ad2b36e6
......@@ -87,11 +87,10 @@ gLevelsetPlane::gLevelsetPlane (const double * pt, const double *norm, int &tag)
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);
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);
}
/*
......@@ -116,7 +115,7 @@ void gLevelsetQuadric::Ax (const double x[3], double res[3], double 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);
+ 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]){
......@@ -168,15 +167,15 @@ void gLevelsetQuadric::computeRotationMatrix (const double dir[3], double t[3][3
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]){
......@@ -186,21 +185,22 @@ void gLevelsetQuadric::computeRotationMatrix (const double dir1[3], const double
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]));
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(){
......@@ -212,10 +212,12 @@ void gLevelsetQuadric::init(){
}
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);
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;
......@@ -225,7 +227,8 @@ gLevelsetGenCylinder :: gLevelsetGenCylinder (const double * pt ,const double *
translate(pt);
}
gLevelsetEllipsoid::gLevelsetEllipsoid (const double *pt, const double *dir, const double &a, const double &b, const double &c, int &tag) : gLevelsetQuadric(tag) {
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);
......@@ -267,7 +270,8 @@ gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *d
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 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));
......@@ -326,15 +330,19 @@ gLevelsetConrod::gLevelsetConrod (const double *pt, const double *dir1, const do
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 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 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 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};
......
contrib/DiscreteIntegration/DILevelset.h
+ 47
− 48
View file @ ad2b36e6
......@@ -202,8 +202,7 @@ 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);
}
......
contrib/DiscreteIntegration/Integration3D.cpp
+ 14
− 11
View file @ ad2b36e6
......@@ -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);
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment