diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index 51a4108dd9f53e40e185fb328fedd67b6d1ebda6..8b0f5fec164d2146122b32e83bb7425b168f0319 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -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 "
@@ -53,9 +54,9 @@ std::string GMSH_DistancePlugin::getHelp() const
"Define the elementary 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"
- "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
@@ -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);
- bbox += SPoint3(v->x(),v->y(),v->z());
+ 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,126 +270,181 @@ 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);
-
- int numNodes = e->getNumVertices();
- std::vector<double> x(numNodes), y(numNodes), z(numNodes);
- std::vector<double> *out = data->incrementList(1, e->getType());
- for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->x());
- for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->y());
- for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->z());
-
- if (maxDim == 3) fprintf(f4,"SS(");
- else if (maxDim == 2) fprintf(f4,"ST(");
- std::vector<double> dist;
- for(int j = 0; j < numNodes; 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*, double>::iterator it = distance_map_PDE.find(v);
- dist.push_back(it->second);
- }
- fprintf(f4,"){");
- for (unsigned int i = 0; i < dist.size(); i++){
- out->push_back(dist[i]);
- if (i) fprintf(f4,",%g", dist[i]);
- else fprintf(f4,"%g", dist[i]);
- }
- fprintf(f4,"};\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> *out = data->incrementList(1, e->getType());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->x());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->y());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->z());
+
+ if (maxDim == 3) fprintf(f4,"SS(");
+ else if (maxDim == 2) fprintf(f4,"ST(");
+ std::vector<double> dist;
+ for(int j = 0; j < numNodes; 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*, double>::iterator it = distance_map.find(v);
+ dist.push_back(it->second);
}
+ fprintf(f4,"){");
+ for (unsigned int i = 0; i < dist.size(); i++){
+ out->push_back(dist[i]);
+ if (i) fprintf(f4,",%g", dist[i]);
+ else fprintf(f4,"%g", dist[i]);
+ }
+ 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)
+#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
-
- //-------------------------------------------------
+#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;
+
+ 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");
- data->Time.push_back(0);
- data->setFileName(fileName.c_str());
- data->finalize();
+ 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;
}
diff --git a/Plugin/PluginManager.cpp b/Plugin/PluginManager.cpp
index b21b41d993752ed3360bccc8429a017986310715..d567ac00739e1c25741d53ef3b6111311ae5f9b9 100644
--- a/Plugin/PluginManager.cpp
+++ b/Plugin/PluginManager.cpp
@@ -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*>
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index 849c949fe3dbd8e6300d9f1870b6cc52bf284567..20cd31efb732844fa8874ab4a968eef651a73710 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -157,7 +157,7 @@ class femTerm {
void addToRightHandSide(dofManager<dataVec> &dm, groupOfElements &C) const
{
- groupOfElements::elementContainer::const_iterator it = C.begin();
+ groupOfElements::elementContainer::const_iterator it = C.begin();
for ( ; it != C.end(); ++it){
MElement *eL = *it;
SElement se(eL);
diff --git a/Solver/groupOfElements.cpp b/Solver/groupOfElements.cpp
index 2ca9fe4c3994d4f62784b269339ff6be99b70824..5da5342ab5ee6f980c352c6f68d4c8c5c9510f5d 100644
--- a/Solver/groupOfElements.cpp
+++ b/Solver/groupOfElements.cpp
@@ -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);
diff --git a/Solver/groupOfElements.h b/Solver/groupOfElements.h
index 8af0d6e0e9811e4c8a512bd629ed74c05d47f8eb..db8f6e0421fdc64d600509a173310897ab768594 100644
--- a/Solver/groupOfElements.h
+++ b/Solver/groupOfElements.h
@@ -37,6 +37,8 @@ class groupOfElements {
}
groupOfElements (GFace*);
+ groupOfElements (GRegion*);
+ groupOfElements(std::vector<MElement*> &elems);
virtual void addPhysical(int dim, int physical) {
elementFilterTrivial filter;
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index 64b6d1890b7764dc9402d7c7b17f9b741f1f0e6a..5946d62e6895ddf24dc427d47ce26ffbe0c76031 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -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);
+
}
};
diff --git a/Solver/orthogonalTerm.h b/Solver/orthogonalTerm.h
index cf1d3c32749c8aff69635e18ff4c679de0516670..6a5bd80ac95bd8125fea1a404b982db2914a50b6 100644
--- a/Solver/orthogonalTerm.h
+++ b/Solver/orthogonalTerm.h
@@ -7,21 +7,29 @@
#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 integrationOrder = 2* (e->getPolynomialOrder() - 1);
int npts;
IntPt *GP;
double jac[3][3];
@@ -29,9 +37,9 @@ 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++){
const double u = GP[i].pt[0];
const double v = GP[i].pt[1];
@@ -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