Skip to content
Snippets Groups Projects
Commit 1138708c authored by Emilie Marchandise's avatar Emilie Marchandise
Browse files

Corrected bug in dgFunctionEvaluator 

Plugin Distrance cleaned
parent 51a358ac
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 34 additions and 244 deletions
Plugin/Distance.cpp
+ 34
− 244
View file @ 1138708c
......@@ -9,14 +9,11 @@
#include "Gmsh.h"
#include "GmshConfig.h"
#include "GModel.h"
#include "MElement.h"
#include "Distance.h"
#include "Context.h"
#include "Numeric.h"
#if defined(HAVE_SOLVER)
#include "simpleFunction.h"
#include "distanceTerm.h"
#include "Context.h"
#include "Numeric.h"
#include "dofManager.h"
#include "linearSystemGMM.h"
#include "linearSystemCSR.h"
......@@ -24,7 +21,7 @@
#include "orthogonalTerm.h"
#include "laplaceTerm.h"
#include "crossConfTerm.h"
#endif
StringXNumber DistanceOptions_Number[] = {
{GMSH_FULLRC, "PhysPoint", NULL, 0.},
......@@ -40,6 +37,7 @@ StringXString DistanceOptions_String[] = {
{GMSH_FULLRC, "Filename", NULL, "distance.pos"}
};
extern "C"
{
GMSH_Plugin *GMSH_RegisterDistancePlugin()
......@@ -61,21 +59,13 @@ std::string GMSH_DistancePlugin::getHelp() const
return "Plugin(Distance) computes distances to physical entities in "
"a mesh.\n\n"
"Define the physical entities to which the distance is computed. "
"If Point=0, Line=0, and Surface=0, then the distance is computed "
"to all the boundaries of the mesh (edges in 2D and faces in 3D).\n\n"
"Define the physical entities to which the distance is computed. If Point=0, Line=0, and Surface=0, then the distance is computed to all the boundaries of the mesh (edges in 2D and faces in 3D).\n\n"
"Computation<0. computes the geometrical euclidian distance "
"(warning: different than the geodesic distance), and Computation=a>0.0 "
"solves a PDE on the mesh with the diffusion constant mu = a*bbox, with "
"bbox being the max size of the bounding box of the mesh (see paper "
"Legrand 2006).\n\n"
"Computation<0. computes the geometrical euclidian distance (warning: different than the geodesic distance), and Computation=a>0.0 solves a PDE on the mesh with the diffusion constant mu = a*bbox, with bbox being the max size of the bounding box of the mesh (see paper Legrand 2006).\n\n"
"Min Scale and max Scale, scale the distance function. If min Scale<0 "
"and max Scale<0, then no scaling is applied to the distance function.\n\n"
"Min Scale and max Scale, scale the distance function. If min Scale<0 and max Scale<0, then no scaling is applied to the distance function.\n\n"
"Plugin(Distance) creates a new distance view and also saves the view "
"in the fileName.pos file.";
"Plugin(Distance) creates a new distance view and also saves the view in the fileName.pos file.";
}
int GMSH_DistancePlugin::getNbOptions() const
......@@ -162,6 +152,7 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
PView *GMSH_DistancePlugin::execute(PView *v)
{
int id_pt = (int) DistanceOptions_Number[0].def;
int id_line = (int) DistanceOptions_Number[1].def;
int id_face = (int) DistanceOptions_Number[2].def;
......@@ -170,9 +161,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
PView *view = new PView();
_data = getDataList(view);
#if defined(HAVE_SOLVER)
#if defined(HAVE_TAUCS)
#ifdef HAVE_TAUCS
linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#else
linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
......@@ -181,20 +170,13 @@ PView *GMSH_DistancePlugin::execute(PView *v)
lsys->setPrec(5.e-8);
#endif
dofManager<double> * dofView = new dofManager<double>(lsys);
#endif
std::vector<GEntity*> _entities;
GModel::current()->getEntities(_entities);
GEntity* ge = _entities[_entities.size()-1];
int integrationPointTetra[2];
if (type==-100){
integrationPointTetra[0]=1;
integrationPointTetra[1]=4;
}
else{
integrationPointTetra[0]=0;
integrationPointTetra[1]=0;
}
int numnodes = 0;
for(unsigned int i = 0; i < _entities.size()-1; i++)
......@@ -230,36 +212,8 @@ PView *GMSH_DistancePlugin::execute(PView *v)
std::vector<int> isInYarn2;
std::vector<SPoint3> closePts2;
if (type==-100){
index.clear();
distancesE.clear();
closePts.clear();
isInYarn.clear();
isInYarn.reserve(totNumNodes);
closePts.reserve(totNumNodes);
distancesE.reserve(totNumNodes);
index.reserve(totNumNodes);
distances2.clear();
distancesE2.clear();
closePts2.clear();
isInYarn2.clear();
distances2.reserve(totNumNodes);
isInYarn2.reserve(totNumNodes);
closePts2.reserve(totNumNodes);
distancesE2.reserve(totNumNodes);
}
for (int i = 0; i < totNumNodes; i++) {
distances.push_back(1.e22);
if (type==-100){
distancesE.push_back(1.e22);
isInYarn.push_back(0);
closePts.push_back(SPoint3(0.,0.,0.));
distances2.push_back(1.e22);
distancesE2.push_back(1.e22);
isInYarn2.push_back(0);
closePts2.push_back(SPoint3(0.,0.,0.));
}
}
int k = 0;
......@@ -270,44 +224,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
MVertex *v = ge->mesh_vertices[j];
pts.push_back(SPoint3(v->x(), v->y(),v->z()));
_distance_map.insert(std::make_pair(v, 0.0));
if (type==-100){
index.push_back(v->getIndex());
_isInYarn_map.insert(std::make_pair(v, 0));
_distanceE_map.insert(std::make_pair(v, 0.0));
}
pt2Vertex[k] = v;
k++;
}
}
if (type==-100){
double jac[3][3];
for (unsigned int i = 0; i < ge->getNumMeshElements(); i++){
MElement *e = ge->getMeshElement(i);
IntPt *ptsi;
int nptsi;
double uvw[4][3];
e->getIntegrationPoints(e->getPolynomialOrder(),&nptsi, &ptsi);
for(int j = 0; j < 4; j++) {
double xyz[3] = {e->getVertex(j)->x(),
e->getVertex(j)->y(),
e->getVertex(j)->z()};
e->xyz2uvw(xyz, uvw[j]);
}
for(int ip = 0; ip < nptsi; ip++){
const double u = ptsi[ip].pt[0];
const double v = ptsi[ip].pt[1];
const double w = ptsi[ip].pt[2];
const double weight = ptsi[ip].weight;
const double detJ = e->getJacobian(u, v, w, jac);
SPoint3 p;
e->pnt(u, v, w, p);
pts.push_back(p);
}
}
}
// Compute geometrical distance to mesh boundaries
//------------------------------------------------------
if (type < 0.0 ){
......@@ -319,12 +240,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
int gDim = g2->dim();
std::vector<int> phys = g2->getPhysicalEntities();
bool computeForEntity = false;
for(unsigned int k = 0; k< phys.size(); k++){
for(int k = 0; k< phys.size(); k++){
int tagp = phys[k];
if (id_pt==0 && id_line==0 && id_face==0 && gDim==_maxDim-1 )
computeForEntity = true;
else if ( (tagp==id_pt && gDim==0)|| (tagp==id_line && gDim==1) ||
(tagp==id_face && gDim==2) )
else if ( (tagp==id_pt && gDim==0)|| (tagp==id_line && gDim==1) || (tagp==id_face && gDim==2) )
computeForEntity = true;
}
if (computeForEntity){
......@@ -339,79 +259,17 @@ PView *GMSH_DistancePlugin::execute(PView *v)
MVertex *v2 = e->getVertex(1);
SPoint3 p1(v1->x(), v1->y(), v1->z());
SPoint3 p2(v2->x(), v2->y(), v2->z());
if((e->getNumVertices() == 2 && order==1) ||
(e->getNumVertices() == 3 && order==2)){
if (type==-100){
// if ( !((p1.x()==p2.x()) & (p1.y()==p2.y()) & (p1.z()==p2.z())) ){
signedDistancesPointsEllipseLine(iDistances, iDistancesE,
iIsInYarn, iClosePts, pts, p1,p2);
// }
}
else{
if((e->getNumVertices() == 2 && order==1) || (e->getNumVertices() == 3 && order==2)){
signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
}
}
else if(e->getNumVertices() == 3 && order==1){
else if((e->getNumVertices() == 3 and order==1) or (e->getNumVertices() == 6 and order==2)){
MVertex *v3 = e->getVertex(2);
SPoint3 p3 (v3->x(),v3->y(),v3->z());
signedDistancesPointsTriangle(iDistances, iClosePts, pts, p1, p2, p3);
}
for (unsigned int kk = 0; kk< pts.size(); kk++) {
if (type==-100){
if( !((p1.x()==p2.x()) & (p1.y()==p2.y()) & (p1.z()==p2.z()))){
if (iIsInYarn[kk]>0){
if (isInYarn[kk]==0){
distances[kk] = iDistances[kk];
distancesE[kk]= iDistancesE[kk];
isInYarn[kk] = iIsInYarn[kk];
closePts[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),
iClosePts[kk].z());
}
else{
if (isInYarn[kk]!=iIsInYarn[kk]){
if (isInYarn2[kk]==0){
distances2[kk] = iDistances[kk];
distancesE2[kk]= iDistancesE[kk];
isInYarn2[kk] = iIsInYarn[kk];
closePts2[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),
iClosePts[kk].z());
}
else{
if (isInYarn2[kk]==iIsInYarn[kk]){
if (iDistancesE[kk] < distancesE2[kk]){
distances2[kk] = iDistances[kk];
distancesE2[kk]= iDistancesE[kk];
closePts2[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),
iClosePts[kk].z());
}
}
}
}
else{
if (iDistancesE[kk] < distancesE[kk]){
distances[kk] = iDistances[kk];
distancesE[kk]= iDistancesE[kk];
closePts[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),
iClosePts[kk].z());
}
}
}
}
else{
if (isInYarn[kk]==0){
if (iDistancesE[kk] < distancesE[kk]){
distances[kk] = iDistances[kk];
distancesE[kk]= iDistancesE[kk];
closePts[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),
iClosePts[kk].z());
}
}
}
}
}
else{
if (fabs(iDistances[kk]) < distances[kk]){
distances[kk] = fabs(iDistances[kk]);
if (std::abs(iDistances[kk]) < distances[kk]){
distances[kk] = std::abs(iDistances[kk]);
MVertex *v = pt2Vertex[kk];
_distance_map[v] = distances[kk];
}
......@@ -419,25 +277,6 @@ PView *GMSH_DistancePlugin::execute(PView *v)
}
}
}
}
if (type==-100){
for (unsigned int kk = 0; kk< pts.size(); kk++) {
if (isInYarn2[kk]>0){
if (distancesE2[kk]>distancesE[kk]){
distances[kk]=distances2[kk];
distancesE[kk]=distancesE2[kk];
isInYarn[kk]=isInYarn2[kk];
closePts[kk]=closePts2[kk];
}
}
if (kk<totNodes){
MVertex *v = pt2Vertex[kk];
_distance_map[v] = distancesE[kk];
_distanceE_map[v] = distances[kk];
_isInYarn_map[v] = isInYarn[kk];
}
}
}
if (!existEntity){
if (id_pt != 0) Msg::Error("The Physical Point does not exist !");
if (id_line != 0) Msg::Error("The Physical Line does not exist !");
......@@ -445,54 +284,6 @@ PView *GMSH_DistancePlugin::execute(PView *v)
return view;
}
printView(_entities, _distance_map);
if (type==-100){
Msg::Info("Writing integrationPointInYarn.pos");
FILE* f5 = fopen("integrationPointInYarn.pos","w");
FILE* f6 = fopen("integrationPointInYarn.bin","wb");
FILE* f7 = fopen("integrationPointInYarn.txt","w");
int j=0;
fprintf(f5,"View \"integrationPointInYarn\"{\n");
for (std::vector<int>::iterator it = isInYarn.begin(); it !=isInYarn.end(); it++) {
if (j>=totNodes){
int iPIY=*it;
fwrite(&iPIY,sizeof(int),1,f6);
fprintf(f7,"%d %lf %lf %lf\n",iPIY,pts[j].x(), pts[j].y(), pts[j].z());
fprintf(f5,"SP(%g,%g,%g){%d};\n",
pts[j].x(), pts[j].y(), pts[j].z(),
*it);
}
j++;
}
fclose(f6);
fclose(f7);
fprintf(f5,"};\n");
fclose(f5);
Msg::Info("Writing isInYarn.pos");
FILE * f4 = fopen("isInYarn.pos","w");
fprintf(f4,"View \"isInYarn\"{\n");
for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
MElement *e = ge->getMeshElement(i);
fprintf(f4,"SS(");
std::vector<int> inYarn;
for(int j = 0; j < e->getNumVertices(); j++) {
MVertex *v = e->getVertex(j);
if(j) fprintf(f4,",%g,%g,%g",v->x(),v->y(), v->z());
else fprintf(f4,"%g,%g,%g", v->x(),v->y(), v->z());
std::map<MVertex*, int>::iterator it = _isInYarn_map.find(v);
inYarn.push_back(it->second);
}
fprintf(f4,"){");
for (unsigned int i=0; i<inYarn.size(); i++){
if (i) fprintf(f4,",%d", inYarn[i]);
else fprintf(f4,"%d", inYarn[i]);
}
fprintf(f4,"};\n");
}
fprintf(f4,"};\n");
fclose(f4);
}
}
// Compute PDE for distance function
......@@ -509,12 +300,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
int gDim = ge->dim();
bool fixForEntity = false;
std::vector<int> phys = ge->getPhysicalEntities();
for(unsigned int k = 0; k< phys.size(); k++){
for(int k = 0; k< phys.size(); k++){
int tagp = phys[k];
if (id_pt==0 && id_line==0 && id_face==0 && gDim==_maxDim-1 )
fixForEntity = true;
else if ( (tagp==id_pt && gDim==0)|| (tagp==id_line && gDim==1) ||
(tagp==id_face && gDim==2) )
else if ( (tagp==id_pt && gDim==0)|| (tagp==id_line && gDim==1) || (tagp==id_face && gDim==2) )
fixForEntity = true;
}
if (fixForEntity){
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment