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

Compute Distance and orthogonal field

parent 023e4a6b
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Plugin/Distance.cpp
+ 192
− 130
View file @ f46002ac
......@@ -18,16 +18,16 @@
#include "linearSystemGMM.h"
#include "linearSystemCSR.h"
#include "linearSystemFull.h"
#include "orthogonalTerm.h"
#include "laplaceTerm.h"
#include "crossConfTerm.h"
#if defined(HAVE_ANN)
#include "ANN/ANN.h"
#endif
StringXNumber DistanceOptions_Number[] = {
{GMSH_FULLRC, "Point", NULL, 0.},
{GMSH_FULLRC, "Line", NULL, 0.},
{GMSH_FULLRC, "Surface", NULL, 0.},
{GMSH_FULLRC, "Computation", NULL, -1.0},
{GMSH_FULLRC, "Computation", NULL, 1.0},
//{GMSH_FULLRC, "Scale", NULL, 0.},
//{GMSH_FULLRC, "Min Scale", NULL, 1.e-3},
//{GMSH_FULLRC, "Max Scale", NULL, 1}
......@@ -46,6 +46,7 @@ extern "C"
}
}
std::string GMSH_DistancePlugin::getHelp() const
{
return "Plugin(Distance) computes distances to elementary entities in "
......@@ -87,26 +88,36 @@ PView *GMSH_DistancePlugin::execute(PView *v)
int id_face = (int) DistanceOptions_Number[2].def;
double type = (double) DistanceOptions_Number[3].def;
PView *view = new PView();
PViewDataList *data = getDataList(view);
#ifdef HAVE_TAUCS
linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#else
linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
lsys->setNoisy(1);
lsys->setGmres(1);
lsys->setPrec(5.e-8);
#endif
dofManager<double> * dofView = new dofManager<double>(lsys);
std::vector<GEntity*> entities;
GModel::current()->getEntities(entities);
int numnodes = 0;
for(unsigned int i = 0; i < entities.size() - 1; i++)
numnodes += entities[i]->mesh_vertices.size();
int totNumNodes = numnodes + entities[entities.size() - 1]->mesh_vertices.size();
int totNumNodes = 0;
for(unsigned int i = 0; i < entities.size() ; i++)
totNumNodes += entities[i]->mesh_vertices.size();
std::map<MVertex*,double* > distance_map_GEOM;
std::map<MVertex*,double > distance_map_PDE;
std::map<MVertex*,double > distance_map;
std::vector<SPoint3> pts;
std::vector<double> distances;
std::vector<MVertex* > pt2Vertex;
pts.clear();
distances.clear();
distances.reserve(totNumNodes);
pt2Vertex.clear();
pts.reserve(totNumNodes);
distances.reserve(totNumNodes);
pt2Vertex.reserve(totNumNodes);
for (int i = 0; i < totNumNodes; i++) distances.push_back(1.e22);
int k = 0;
......@@ -117,8 +128,8 @@ PView *GMSH_DistancePlugin::execute(PView *v)
for(unsigned int j = 0; j < ge->mesh_vertices.size(); j++){
MVertex *v = ge->mesh_vertices[j];
pts.push_back(SPoint3(v->x(), v->y(),v->z()));
distance_map_GEOM.insert(std::make_pair(v, &(distances[k])));
distance_map_PDE.insert(std::make_pair(v, 0.0));
distance_map.insert(std::make_pair(v, 0.0));
pt2Vertex[k] = v;
k++;
}
}
......@@ -154,8 +165,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
signedDistancesPointsTriangle(iDistances, iClosePts, pts, p1, p2, p3);
}
for (unsigned int kk = 0; kk< pts.size(); kk++) {
if (std::abs(iDistances[kk]) < distances[kk])
if (std::abs(iDistances[kk]) < distances[kk]){
distances[kk] = std::abs(iDistances[kk]);
MVertex *v = pt2Vertex[kk];
distance_map[v] = distances[kk];
}
}
}
}
......@@ -170,7 +184,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
// for(unsigned int i = 0; i < g2->getNumMeshElements(); i++)
data->setName("distance GEOM");
Msg::Info("Writing %g", fileName.c_str());
Msg::Info("Writing %s", fileName.c_str());
FILE * f3 = fopen( fileName.c_str(),"w");
fprintf(f3,"View \"distance GEOM\"{\n");
for(unsigned int ii = 0; ii < entities.size(); ii++){
......@@ -192,8 +206,8 @@ PView *GMSH_DistancePlugin::execute(PView *v)
MVertex *v = e->getVertex(j);
if(j) fprintf(f3,",%g,%g,%g",v->x(),v->y(), v->z());
else fprintf(f3,"%g,%g,%g", v->x(),v->y(), v->z());
std::map<MVertex*, double*>::iterator it = distance_map_GEOM.find(v);
dist.push_back(*(it->second));
std::map<MVertex*, double>::iterator it = distance_map.find(v);
dist.push_back(it->second);
}
fprintf(f3,"){");
for (unsigned int i = 0; i < dist.size(); i++){
......@@ -215,17 +229,6 @@ PView *GMSH_DistancePlugin::execute(PView *v)
#if defined(HAVE_SOLVER)
#ifdef HAVE_TAUCS
linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#else
linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
lsys->setNoisy(1);
lsys->setGmres(1);
lsys->setPrec(5.e-8);
#endif
dofManager<double> myAssembler(lsys);
SBoundingBox3d bbox;
for(unsigned int i = 0; i < entities.size(); i++){
GEntity* ge = entities[i];
......@@ -237,21 +240,25 @@ PView *GMSH_DistancePlugin::execute(PView *v)
else if ( (tag==id_pt && gDim==0)|| (tag==id_line && gDim==1) || (tag==id_face && gDim==2) )
fixForEntity = true;
if (fixForEntity){
for(unsigned int j = 0; j < ge->mesh_vertices.size(); j++){
MVertex *v = ge->mesh_vertices[j];
myAssembler.fixVertex(v, 0, 1, 0.0);
for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
MElement *t = ge->getMeshElement(i);
for(int k = 0; k < t->getNumVertices(); k++){
MVertex *v = t->getVertex(k);
dofView->fixVertex(v, 0, 1, 0.0);
bbox += SPoint3(v->x(),v->y(),v->z());
}
}
}
}
std::vector<MElement *> allElems;
for(unsigned int ii = 0; ii < entities.size(); ii++){
if(entities[ii]->dim() == maxDim) {
GEntity *ge = entities[ii];
for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
MElement *t = ge->getMeshElement(i);
allElems.push_back(t);
for(int k = 0; k < t->getNumVertices(); k++){
myAssembler.numberVertex(t->getVertex(k), 0, 1);
dofView->numberVertex(t->getVertex(k), 0, 1);
}
}
}
......@@ -263,42 +270,33 @@ PView *GMSH_DistancePlugin::execute(PView *v)
simpleFunction<double> DIFF(mu*mu), ONE(1.0);
distanceTerm distance(GModel::current(), 1, &DIFF, &ONE);
for(unsigned int ii = 0; ii < entities.size(); ii++){
if(entities[ii]->dim() == maxDim) {
GEntity *ge = entities[ii];
for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
SElement se(ge->getMeshElement(i));
distance.addToMatrix(myAssembler, &se);
}
groupOfElements g((GFace*)ge);
distance.addToRightHandSide(myAssembler, g);
}
for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
SElement se((*it));
distance.addToMatrix(*dofView, &se);
}
groupOfElements gr(allElems);
distance.addToRightHandSide(*dofView, gr);
Msg::Info("Distance Computation: Assembly done");
lsys->systemSolve();
Msg::Info("Distance Computation: System solved");
for(std::map<MVertex*,double >::iterator itv = distance_map_PDE.begin();
itv !=distance_map_PDE.end() ; ++itv){
for(std::map<MVertex*,double >::iterator itv = distance_map.begin();
itv !=distance_map.end() ; ++itv){
MVertex *v = itv->first;
double value;
myAssembler.getDofValue(v, 0, 1, value);
dofView->getDofValue(v, 0, 1, value);
value = std::min(0.9999, value);
double dist = -mu * log(1. - value);
itv->second = dist;
}
lsys->clear();
data->setName("distance PDE");
Msg::Info("Writing %d", fileName.c_str());
Msg::Info("Writing %s", fileName.c_str());
FILE * f4 = fopen(fileName.c_str(),"w");
fprintf(f4,"View \"distance PDE\"{\n");
for(unsigned int ii = 0; ii < entities.size(); ii++){
if(entities[ii]->dim() == maxDim) {
for(unsigned int i = 0; i < entities[ii]->getNumMeshElements(); i++){
MElement *e = entities[ii]->getMeshElement(i);
for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
MElement *e = *it;
int numNodes = e->getNumVertices();
std::vector<double> x(numNodes), y(numNodes), z(numNodes);
std::vector<double> *out = data->incrementList(1, e->getType());
......@@ -313,7 +311,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
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*, double>::iterator it = distance_map_PDE.find(v);
std::map<MVertex*, double>::iterator it = distance_map.find(v);
dist.push_back(it->second);
}
fprintf(f4,"){");
......@@ -324,65 +322,129 @@ PView *GMSH_DistancePlugin::execute(PView *v)
}
fprintf(f4,"};\n");
}
}
}
fprintf(f4,"};\n");
fclose(f4);
#endif
}
data->Time.push_back(0);
data->setFileName(fileName.c_str());
data->finalize();
//compute also orthogonal vector to distance field
// A Uortho = -C DIST
//------------------------------------------------
// #if defined(HAVE_SOLVER)
// #ifdef HAVE_TAUCS
// linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
// #else
// linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
// lsys->setNoisy(1);
// lsys->setGmres(1);
// lsys->setPrec(5.e-8);
// #endif
// dofManager<double> myAssembler(lsys);
// for(unsigned int ii = 0; ii < entities.size(); ii++){
// if(entities[ii]->dim() == maxDim) {
// GEntity *ge = entities[ii];
// for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
// MElement *t = ge->getMeshElement(i);
// for(int k = 0; k < t->getNumVertices(); k++){
// myAssembler.numberVertex(t->getVertex(k), 0, 1);
// }
// }
// }
// }
// simpleFunction<double> ONE(1.0);
// simpleFunction<double> MONE(-1.0 );
// laplaceTerm laplace(model(), 1, &ONE);
// crossConfTerm cross12(model(), 1, 1, &ONE);
// for(unsigned int ii = 0; ii < entities.size(); ii++){
// if(entities[ii]->dim() == maxDim) {
// GEntity *ge = entities[ii];
// for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
// SElement se(ge->getMeshElement(i));
// laplace.addToMatrix(myAssembler, &se);
// }
// //groupOfElements g((GFace*)ge); //WARNING GFACE HE
// //distance.addToRightHandSide(myAssembler, g);
// }
// }
// #endif
#if defined(HAVE_SOLVER)
//-------------------------------------------------
#ifdef HAVE_TAUCS
linearSystemCSRTaucs<double> *lsys2 = new linearSystemCSRTaucs<double>;
#else
linearSystemCSRGmm<double> *lsys2 = new linearSystemCSRGmm<double>;
lsys->setNoisy(1);
lsys->setGmres(1);
lsys->setPrec(5.e-8);
#endif
dofManager<double> myAssembler(lsys2);
simpleFunction<double> ONE(1.0);
double dMax = 0.03;
data->Time.push_back(0);
data->setFileName(fileName.c_str());
data->finalize();
std::vector<MElement *> allElems;
for(unsigned int ii = 0; ii < entities.size(); ii++){
if(entities[ii]->dim() == maxDim) {
GEntity *ge = entities[ii];
for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
MElement *t = ge->getMeshElement(i);
double vMean = 0.0;
for(int k = 0; k < t->getNumVertices(); k++){
std::map<MVertex*, double>::iterator it = distance_map.find(t->getVertex(k));
vMean += it->second;
}
vMean/=t->getNumVertices();
if (vMean < dMax)
allElems.push_back(ge->getMeshElement(i));
}
}
}
int mid = (int)floor(allElems.size()/2);
MElement *e = allElems[mid];
MVertex *vFIX = e->getVertex(0);
myAssembler.fixVertex(vFIX, 0, 1, 0.0);
for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
MElement *t = *it;
for(int k = 0; k < t->getNumVertices(); k++)
myAssembler.numberVertex(t->getVertex(k), 0, 1);
}
orthogonalTerm *ortho;
ortho = new orthogonalTerm(GModel::current(), 1, &ONE, &distance_map);
// if (type < 0)
// ortho = new orthogonalTerm(GModel::current(), 1, &ONE, view);
// else
// ortho = new orthogonalTerm(GModel::current(), 1, &ONE, dofView);
for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
SElement se((*it));
ortho->addToMatrix(myAssembler, &se);
}
groupOfElements gr(allElems);
ortho->addToRightHandSide(myAssembler, gr);
Msg::Info("Orthogonal Computation: Assembly done");
lsys2->systemSolve();
Msg::Info("Orthogonal Computation: System solved");
PView *view2 = new PView();
PViewDataList *data2 = getDataList(view2);
data2->setName("ortogonal field");
Msg::Info("Writing orthogonal.pos");
FILE * f5 = fopen("orthogonal.pos","w");
fprintf(f5,"View \"orthogonal\"{\n");
for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
MElement *e = *it;
int numNodes = e->getNumVertices();
std::vector<double> x(numNodes), y(numNodes), z(numNodes);
std::vector<double> *out2 = data2->incrementList(1, e->getType());
for(int nod = 0; nod < numNodes; nod++) out2->push_back((e->getVertex(nod))->x());
for(int nod = 0; nod < numNodes; nod++) out2->push_back((e->getVertex(nod))->y());
for(int nod = 0; nod < numNodes; nod++) out2->push_back((e->getVertex(nod))->z());
if (maxDim == 3) fprintf(f5,"SS(");
else if (maxDim == 2) fprintf(f5,"ST(");
std::vector<double> orth;
for(int j = 0; j < numNodes; j++) {
MVertex *v = e->getVertex(j);
if(j) fprintf(f5,",%g,%g,%g",v->x(),v->y(), v->z());
else fprintf(f5,"%g,%g,%g", v->x(),v->y(), v->z());
double value;
myAssembler.getDofValue(v, 0, 1, value );
orth.push_back(value);
}
fprintf(f5,"){");
for (unsigned int i = 0; i < orth.size(); i++){
out2->push_back(orth[i]);
if (i) fprintf(f5,",%g", orth[i]);
else fprintf(f5,"%g", orth[i]);
}
fprintf(f5,"};\n");
}
fprintf(f5,"};\n");
fclose(f5);
lsys->clear();
lsys2->clear();
data2->Time.push_back(0);
data2->setFileName("orthogonal.pos");
data2->finalize();
#endif
//-------------------------------------------------
return view;
}
Plugin/PluginManager.cpp
+ 3
− 0
View file @ f46002ac
......@@ -29,6 +29,7 @@
#include "Curl.h"
#include "Divergence.h"
#include "Annotate.h"
#include "Distance.h"
#include "Remove.h"
#include "MakeSimplex.h"
#include "Smooth.h"
......@@ -198,6 +199,8 @@ void PluginManager::registerDefaultPlugins()
("Curl", GMSH_RegisterCurlPlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("Divergence", GMSH_RegisterDivergencePlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("Distance", GMSH_RegisterDistancePlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("Annotate", GMSH_RegisterAnnotatePlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
......
Solver/femTerm.h
+ 1
− 1
View file @ f46002ac
Solver/groupOfElements.cpp
+ 18
− 0
View file @ f46002ac
......@@ -8,6 +8,24 @@ groupOfElements::groupOfElements(GFace*gf)
addElementary(gf, filter);
}
groupOfElements::groupOfElements(GRegion*gr)
{
elementFilterTrivial filter;
addElementary(gr, filter);
}
groupOfElements::groupOfElements(std::vector<MElement*> &elems)
{
elementFilterTrivial filter;
for (std::vector<MElement*>::iterator it = elems.begin(); it != elems.end(); it++){
MElement *e = *it;
if (filter(e)){
insert(e);
}
}
}
void groupOfElements::addElementary(GEntity *ge, const elementFilter &filter){
for (unsigned int j = 0; j < ge->getNumMeshElements(); j++){
MElement *e = ge->getMeshElement(j);
......
Solver/groupOfElements.h
+ 2
− 0
View file @ f46002ac
......@@ -37,6 +37,8 @@ class groupOfElements {
}
groupOfElements (GFace*);
groupOfElements (GRegion*);
groupOfElements(std::vector<MElement*> &elems);
virtual void addPhysical(int dim, int physical) {
elementFilterTrivial filter;
......
Solver/helmholtzTerm.h
+ 2
− 0
View file @ f46002ac
......@@ -47,6 +47,7 @@ class helmholtzTerm : public femTerm<scalar> {
}
virtual void elementMatrix(SElement *se, fullMatrix<scalar> &m) const
{
MElement *e = se->getMeshElement();
// compute integration rule
const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() :
......@@ -98,6 +99,7 @@ class helmholtzTerm : public femTerm<scalar> {
for (int j = 0; j < nbNodes; j++)
for (int k = 0; k < j; k++)
m(k, j) = m(j, k);
}
};
......
Solver/orthogonalTerm.h
+ 44
− 17
View file @ f46002ac
......@@ -7,20 +7,28 @@
#define _ORTHOGONAL_TERM_H_
#include "helmholtzTerm.h"
#include "dofManager.h"
class orthogonalTerm : public helmholtzTerm<double> {
protected:
fullVector<double> *_value;
PView *_pview;
dofManager<double> *_dofView;
bool withDof;
std::map<MVertex*,double > *_distance_map;
public:
orthogonalTerm(GModel *gm, int iField, fullVector<double> &value)
: helmholtzTerm<double>(gm, iField, iField, 1.0, 0), _value(value) {}
orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k, std::map<MVertex*,double > *distance_map)
: helmholtzTerm<double>(gm, iField, iField, k, 0), _distance_map(distance_map), withDof(false) {}
orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k, PView *pview)
: helmholtzTerm<double>(gm, iField, iField, k, 0), _pview(pview), withDof(false) {}
orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k, dofManager<double> *dofView)
: helmholtzTerm<double>(gm, iField, iField, k, 0), _dofView(dofView), withDof(true) {}
void elementVector(SElement *se, fullVector<double> &m) const
{
MElement *e = se->getMeshElement();
int nbNodes = e->getNumVertices();
//fill elementary matrix mat(i,j)
int nbNodes = e->getNumVertices();
int integrationOrder = 2* (e->getPolynomialOrder() - 1);
int npts;
IntPt *GP;
......@@ -29,7 +37,7 @@ class orthogonalTerm : public helmholtzTerm<double> {
SVector3 Grads [256];
double grads[256][3];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
fullMatrix<double> mat;
fullMatrix<double> mat(nbNodes,nbNodes);
mat.setAll(0.);
for (int i = 0; i < npts; i++){
......@@ -39,7 +47,6 @@ class orthogonalTerm : public helmholtzTerm<double> {
const double weight = GP[i].weight;
const double detJ = e->getJacobian(u, v, w, jac);
SPoint3 p; e->pnt(u, v, w, p);
const double _diff = (*_diffusivity)(p.x(), p.y(), p.z());
inv3x3(jac, invjac);
e->getGradShapeFunctions(u, v, w, grads);
for (int j = 0; j < nbNodes; j++){
......@@ -51,25 +58,45 @@ class orthogonalTerm : public helmholtzTerm<double> {
invjac[2][2] * grads[j][2]);
}
SVector3 N (jac[2][0], jac[2][1], jac[2][2]);
for (int j = 0; j < nbNodes; j++){
for (int k = 0; k <= j; k++){
mat(j, k) += dot(crossprod(Grads[j], Grads[k]), N) * weight * detJ * _diff;
}
}
for (int j = 0; j < nbNodes; j++)
for (int k = 0; k <= j; k++)
mat(j, k) += dot(crossprod(Grads[j], Grads[k]), N) * weight * detJ;
}
for (int j = 0; j < nbNodes; j++)
for (int k = 0; k < j; k++)
mat(k, j) = -1.* m(j, k);
mat(k, j) = -1.* mat(j, k);
//2) compute vector m(i) = mat(i,j)*val(j)
fullVector<double> val(nbNodes);
val.scale(0.);
for (int i = 0; i < nbNodes; i++){
std::map<MVertex*, double>::iterator it = _distance_map->find(e->getVertex(i));
val(i) = it->second;
}
/* if( withDof){ */
/* for (int i = 0; i < nbNodes; i++){ */
/* _dofView->getDofValue( e->getVertex(i), 0, 1, val(i)); */
/* printf("val=%g \n", val(i)); */
/* } */
/* } */
/* else{ */
/* printf("with no dof \n"); */
/* exit(1); */
/* /\* PViewData *data = view->getData(); *\/ */
/* /\* for (int i = 0; i < nbNodes; i++){ *\/ */
/* /\* data->getValue(0, e->, e->getNum(), e->getVertex(i)->getNum(), nod, 0, val(i)); *\/ */
/* /\* } *\/ */
/* } */
m.scale(0.);
for (int i = 0; i < nbNodes; i++)
for (int j = 0; j < nbNodes; j++)
m(i) += mat(i,j)*val(j);
m(i) += -mat(i,j)*val(j);
}
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment