diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index b32a78c27f230451057e8dff1f1babb4a2776542..a5f4755ebf7b713fb05588de25890f465634312c 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2278,7 +2278,7 @@ void GModel::createTopologyFromMesh(int ignoreHoles)
makeDiscreteFacesSimplyConnected();
// // TEST !!!!!!!!
- if (0)
+ if (1)
{
createTopologyFromMeshNew ();
double t2 = Cpu();
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 8f493f8cf2fbdb6aa59a9e1656e4fc55bdf6493c..a4dfab09db998b2d609d818c99d753ced3bc0410 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -27,6 +27,11 @@
#include "qualityMeasuresJacobian.h"
#endif
+#if defined(HAVE_MUMPS)
+#include "linearSystemMUMPS.h"
+#endif
+
+
#if defined(HAVE_OPTHOM)
#include "OptHomRun.h"
#include "MeshQualityOptimizer.h"
@@ -182,7 +187,7 @@ static int discreteDiskFaceInEle(void *a, double*c)// # mark
double Xi[2];
double U[2] = {c[0],c[1]};
bool ok = uv2xi(t,U,Xi);
- double eps = 1e-8;
+ double eps = 1e-6;
if(ok && Xi[0] > -eps && Xi[1] > -eps && 1. - Xi[0] - Xi[1] > -eps)
return 1;
@@ -284,21 +289,10 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
_totLength = geoTriangulation->bord.rbegin()->first;
_U0 = geoTriangulation->bord.rbegin()->second;
orderVertices(_totLength, _U0, _coords);
- parametrize(false);
+ parametrize();
buildOct(CAD);
printParamMesh();
- if (!checkOrientationUV()){
- Msg::Info("discreteDiskFace:: parametrization is not one-to-one; it is going to be corrected.");
- if(_order==1)
- parametrize(true);
- else
- optimize();
- buildOct(CAD);
- printParamMesh();
- }
- putOnView(true,false);
- if(!checkOrientationUV()) Msg::Fatal("discreteDiskFace:: failing to compute a one-to-one mapping");
- else Msg::Info("Parameterization done :-)");
+ Msg::Info("Parameterization done (GFace %d)",gf->tag());
}
/*end BUILDER*/
@@ -339,31 +333,73 @@ void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
Octree_Arrange(oct);
}
-bool discreteDiskFace::parametrize(bool one2one)// const
-{ // #improveme
+/*
+void discreteDiskFace::OneToOneParametrization () {
+ v2t_cont adj;
+ buildVertexToTriangle(discrete_triangles, adj);
+ int count=0;
+ {
+ v2t_cont :: iterator it = adj.begin();
+ while (it != adj.end()) {
+ it->first->setIndex(count++);
+ ++it;
+ }
+ }
+ double *rhs = new double [2*count];
+ double *x = new double [2*count];
+ {
+ for(size_t i = 0; i < _U0.size(); i++){
+ MVertex *v = _U0[i];
+ const double theta = 2 * M_PI * _coords[i];
+ myAssemblerU.fixVertex(v, 0, 1,cos(theta));
+ myAssemblerV.fixVertex(v, 0, 1,sin(theta));
+ }
+ }
+ {
+ v2t_cont :: iterator it = adj.begin();
+ while (it != adj.end()) {
+ it->first->setIndex(i++);
+ ++it;
+ }
+ }
+
+ std::vector<int> ij;
+ std::vector<double> val;
+
+}
+*/
- linearSystemCSRTaucs<double> *lsys_u;
- linearSystemCSRTaucs<double> *lsys_v;
+bool discreteDiskFace::parametrize()// const
+{ // #improveme
- lsys_u = new linearSystemCSRTaucs<double>; // taucs :-)
- lsys_v = new linearSystemCSRTaucs<double>; // taucs :-)
+ linearSystem<double> * lsys_u, *lsys_v;
+
+#ifdef HAVE_MUMPS
+ lsys_u = new linearSystemMUMPS<double>;
+ lsys_v = new linearSystemMUMPS<double>;
+#else
+ linearSystemCSRGmm<double> * lsys_u1 = new linearSystemCSRGmm<double>;
+ linearSystemCSRGmm<double> * lsys_v1 = new linearSystemCSRGmm<double>;
+ lsys_u1->setGmres(1);
+ lsys_v1->setGmres(1);
+ lsys_u=lsys_u1;
+ lsys_v=lsys_v1;
+#endif
dofManager<double> myAssemblerU(lsys_u); // hashing
dofManager<double> myAssemblerV(lsys_v);
-
- if(!one2one){
- for(size_t i = 0; i < _U0.size(); i++){
- MVertex *v = _U0[i];
- const double theta = 2 * M_PI * _coords[i];
- if(i%_order==0){
- myAssemblerU.fixVertex(v, 0, 1,cos(theta));
- myAssemblerV.fixVertex(v, 0, 1,sin(theta));
- }
- else{//#TEST
- myAssemblerU.fixVertex(v, 0, 1,1./_order*((_order-(i%_order))*cos(2*M_PI*_coords[i-(i%_order)])+(i%_order)*cos(2*M_PI*_coords[i+(_order-(i%_order))])));
- myAssemblerV.fixVertex(v, 0, 1,1./_order*((_order-(i%_order))*sin(2*M_PI*_coords[i-(i%_order)])+(i%_order)*sin(2*M_PI*_coords[i+(_order-(i%_order))])));
- }
+
+ for(size_t i = 0; i < _U0.size(); i++){
+ MVertex *v = _U0[i];
+ const double theta = 2 * M_PI * _coords[i];
+ if(i%_order==0){
+ myAssemblerU.fixVertex(v, 0, 1,cos(theta));
+ myAssemblerV.fixVertex(v, 0, 1,sin(theta));
+ }
+ else{//#TEST
+ myAssemblerU.fixVertex(v, 0, 1,1./_order*((_order-(i%_order))*cos(2*M_PI*_coords[i-(i%_order)])+(i%_order)*cos(2*M_PI*_coords[i+(_order-(i%_order))])));
+ myAssemblerV.fixVertex(v, 0, 1,1./_order*((_order-(i%_order))*sin(2*M_PI*_coords[i-(i%_order)])+(i%_order)*sin(2*M_PI*_coords[i+(_order-(i%_order))])));
}
}
@@ -382,8 +418,8 @@ bool discreteDiskFace::parametrize(bool one2one)// const
double t1 = Cpu();
simpleFunction<double> ONE(1.0);
- laplaceTerm mappingU(0, 1, &ONE);
- laplaceTerm mappingV(0, 1, &ONE);
+ convexCombinationTerm mappingU(0, 1, &ONE);
+ convexCombinationTerm mappingV(0, 1, &ONE);
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
SElement se(discrete_triangles[i]);
@@ -391,13 +427,6 @@ bool discreteDiskFace::parametrize(bool one2one)// const
mappingV.addToMatrix(myAssemblerV, &se);
}
-
- if(one2one){
- Msg::Info("discreteDiskFace::parametrize \t Modifying discrete system");
- checklsys(lsys_u,&myAssemblerU,1);
- checklsys(lsys_v,&myAssemblerV,0);
- }
-
double t2 = Cpu();
Msg::Debug("Assembly done in %8.3f seconds", t2 - t1);
lsys_u->systemSolve();
@@ -429,97 +458,30 @@ bool discreteDiskFace::parametrize(bool one2one)// const
return true;
}
-void discreteDiskFace::checklsys(linearSystemCSRTaucs<double>* lsys,dofManager<double>* myAssembler,int U)
-{
-
- INDEX_TYPE *ai;
- INDEX_TYPE *jptr;
- double *a;
- //double b;
- lsys->getMatrix(jptr,ai,a);
- int NbUnk = lsys->getNbUnk();
- // int NNZ = lsys->getNNZ();
- double tr = 0.;
- //printf("------------------------%d\n######################%d\n______________________\n",NNZ,NbUnk);
- for(int i=0; i<NbUnk; i++){
- tr += a[jptr[i]];
- /*lsys->getFromRightHandSide(i,b);
- printf("%d: %d->%d\n",i,jptr[i],jptr[i+1]);
- for(int j=jptr[i]; j<jptr[i+1]; j++)
- printf("%d: %f \t",ai[j],a[j]);
- printf("=%f \n",b);*/
- }
- tr /= NbUnk;
-
- double alpha = .1;
- for(int i=0; i<NbUnk; i++)
- for(int j=jptr[i]+1; j<jptr[i+1]; j++)
- if(a[j] > -alpha*tr){
- double delta = a[j] + alpha*tr;
- a[jptr[i]] += delta;
- a[jptr[ai[j]]] += delta;
- a[j] -= delta;
- }
- /*
- printf("\n");
- for(int i=0; i<NbUnk; i++){
- lsys->getFromRightHandSide(i,b);
- printf("%d: %d->%d\n",i,jptr[i],jptr[i+1]);
- for(int j=jptr[i]; j<jptr[i+1]; j++)
- printf("%d: %f \t",ai[j],a[j]);
- printf("=%f \n",b);
- }
- */
- std::map<int,int> dirichlet;
- for(unsigned int i=0; i<_U0.size(); i++)
- dirichlet[myAssembler->getDofNumber(Dof(_U0[i]->getNum(),10000))] = i;
-
-
- for(int i=0; i<NbUnk; i++){
- if(dirichlet.find(i)!=dirichlet.end()){
- lsys->addToMatrix(i,i,-a[jptr[i]]+1);
- double theta = 2 * M_PI * _coords[dirichlet[i]];
- double val = U*cos(theta) + (1-U)*sin(theta);
- lsys->addToRightHandSide(i,val);
- for(int j=jptr[i]+1; j<jptr[i+1]; j++){
- if(dirichlet.find(ai[j])==dirichlet.end())
- lsys->addToRightHandSide(ai[j],-a[j]*val);
- a[j] = 0;
- }
- }
- else{
- for(int j=jptr[i]+1; j<jptr[i+1]; j++){
- if(dirichlet.find(ai[j])!=dirichlet.end()){
- double theta = 2 * M_PI * _coords[dirichlet[ai[j]]];
- double val = U*cos(theta) + (1-U)*sin(theta);
- lsys->addToRightHandSide(i,-a[j]*val);
- a[j] = 0;
- }
- }
- }//end else
- }//end for i
- /*
- lsys->getMatrix(jptr,ai,a);
- NbUnk = lsys->getNbUnk();
- NNZ = lsys->getNNZ();
- printf("--#--#--#--\n");
- for(int i=0; i<NbUnk; i++){
- lsys->getFromRightHandSide(i,b);
- printf("%d: %d->%d\n",i,jptr[i],jptr[i+1]);
- for(int j=jptr[i]; j<jptr[i+1]; j++)
- printf("%d: %f \t",ai[j],a[j]);
- printf("=%f \n",b);
- }
- */
-
-}
+//void discreteDiskFace::checklsys(linearSystemCSR<double>* lsys,dofManager<double>* myAssembler,int U)
+//{
+//}
void discreteDiskFace::getTriangleUV(const double u,const double v,
discreteDiskFaceTriangle **mt,
double &_xi, double &_eta)const{
double uv[3] = {u,v,0.};
*mt = (discreteDiskFaceTriangle*) Octree_Search(uv,oct);
- if (!(*mt)) return;
+ if (!(*mt)) {
+ for (unsigned int i=0; i<discrete_triangles.size()-geoTriangulation->fillingHoles.size(); i++){
+ discreteDiskFaceTriangle *ct = &_ddft[i];
+ double Xi[2];
+ int xxx = discreteDiskFaceInEle(ct, Xi);
+ if (xxx) {
+ *mt = ct;
+ _xi = Xi[0];
+ _eta = Xi[1];
+ return;
+ }
+ }
+ Msg::Debug("discreteDiskFace::getTriangleUV(), didn't find the reference coordinate (xi;eta) for (u;v)=(%f;%f) among %d triangles",u,v,discrete_triangles.size()-geoTriangulation->fillingHoles.size());
+ return;
+ }
double Xi[2];
double U[2] = {u,v};
@@ -841,7 +803,7 @@ void discreteDiskFace::secondDer(const SPoint2 ¶m,
return;
}
-void discreteDiskFace::putOnView(bool Xu, bool Ux)
+void discreteDiskFace::putOnView(int iFace, int iMap, bool Xu, bool Ux)
{// #improveme using built-in methods
char mybuffer [64];
@@ -851,28 +813,27 @@ void discreteDiskFace::putOnView(bool Xu, bool Ux)
if(Xu){
-
- sprintf(mybuffer, "param_u_part%d_order%d.pos",
- initialTriangulation->idNum,_order);
+ sprintf(mybuffer, "param_u_gface%d_part%d_order%d.pos",
+ iFace, iMap,_order);
view_u = Fopen(mybuffer,"w");
- sprintf(mybuffer, "param_v_part%d_order%d.pos",
- initialTriangulation->idNum,_order);
+ sprintf(mybuffer, "param_v_gface%d_part%d_order%d.pos",
+ iFace, iMap,_order);
view_v = Fopen(mybuffer,"w");
}
if(Ux){
- sprintf(mybuffer, "UVx_part%d_order%d.pos",
- initialTriangulation->idNum,_order);
+ sprintf(mybuffer, "UVx_gface%d_part%d_order%d.pos",
+ iFace, iMap,_order);
UVx = Fopen(mybuffer,"w");
- sprintf(mybuffer, "UVy_part%d_order%d.pos",
- initialTriangulation->idNum,_order);
+ sprintf(mybuffer, "UVy_gface%d_part%d_order%d.pos",
+ iFace, iMap,_order);
UVy = Fopen(mybuffer,"w");
- sprintf(mybuffer, "UVz_part%d_order%d.pos",
- initialTriangulation->idNum,_order);
+ sprintf(mybuffer, "UVz_gface%d_part%d_order%d.pos",
+ iFace, iMap,_order);
UVz = Fopen(mybuffer,"w");
}
@@ -978,10 +939,132 @@ void discreteDiskFace::putOnView(bool Xu, bool Ux)
// useful for mesh generators
// Intersection of a circle and a plane
+//FILE *allProblems = NULL;
+//void openProblemsON(void){
+// allProblems = fopen ("op.pos","w");
+//}
+//void openProblemsOFF(void){
+// fclose(allProblems);
+// allProblems = NULL;
+//}
+
GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
+ const SVector3 &p, const double &R,
+ double uv[2]) const
+{
+ SVector3 n = crossprod(n1,n2);
+ n.normalize();
+ // printf("n %g %g %g\n",n.x(), n.y(), n.z());
+ const int N = (int)(discrete_triangles.size()-geoTriangulation->fillingHoles.size());
+ for (int i=-1;i<N;i++){
+ discreteDiskFaceTriangle *ct = NULL;
+ double U,V;
+ if (i == -1) getTriangleUV(uv[0],uv[1], &ct, U,V);
+ else ct = &_ddft[i];
+ if (!ct) continue;
+ SVector3 v0(ct->tri->getVertex(0)->x(),ct->tri->getVertex(0)->y(),
+ ct->tri->getVertex(0)->z());
+ SVector3 v1(ct->tri->getVertex(1)->x(),ct->tri->getVertex(1)->y(),
+ ct->tri->getVertex(1)->z());
+ SVector3 v2(ct->tri->getVertex(2)->x(),ct->tri->getVertex(2)->y(),
+ ct->tri->getVertex(2)->z());
+ SVector3 t1 = v1 - v0;
+ SVector3 t2 = v2 - v0;
+ SVector3 t = crossprod(t1,t2);
+ t.normalize();
+ SVector3 d = crossprod(n,t);
+ if (d.norm() < 1.e-12) continue;
+ d.normalize();
+ double rhs[2] = {dot(n,p), dot(v0,t)};
+ double r[2];
+ double m[2][2];
+ SVector3 x0(0,0,0);
+ m[0][0] = n.y();
+ m[0][1] = n.z();
+ m[1][0] = t.y();
+ m[1][1] = t.z();
+ if (fabs(det2x2(m)) > 1.e-8){
+ sys2x2(m,rhs,r);
+ x0 = SVector3(0,r[0],r[1]);
+ }
+ else {
+ m[0][0] = n.x();
+ m[0][1] = n.z();
+ m[1][0] = t.x();
+ m[1][1] = t.z();
+ if (fabs(det2x2(m)) > 1.e-8){
+ sys2x2(m,rhs,r);
+ x0 = SVector3(r[0],0,r[1]);
+ }
+ else {
+ m[0][0] = n.x();
+ m[0][1] = n.y();
+ m[1][0] = t.x();
+ m[1][1] = t.y();
+ if (sys2x2(m,rhs,r)) {
+ x0 = SVector3(r[0],r[1],0);
+ }
+ else{
+ printf("mauvaise pioche\n");
+ continue;
+ }
+ }
+ }
+
+ const double a = 1.0;
+ const double b = -2*dot(d,p-x0);
+ const double c = dot(p-x0,p-x0) - R*R;
+ const double delta = b*b-4*a*c;
+ if (delta >= 0){
+ double sign = (dot(n2,d) > 0)? 1.0:-1.0;
+ const double ta = (-b + sign*sqrt(delta)) / (2.*a);
+ const double tb = (-b - sign*sqrt(delta)) / (2.*a);
+ SVector3 s[2] = {x0 + d * ta, x0 + d * tb};
+ for (int IT=0;IT<2;IT++){
+ double mat[2][2], b[2],uv[2];
+ mat[0][0] = dot(t1,t1);
+ mat[1][1] = dot(t2,t2);
+ mat[0][1] = mat[1][0] = dot(t1,t2);
+ b[0] = dot(s[IT]-v0,t1) ;
+ b[1] = dot(s[IT]-v0,t2) ;
+ sys2x2(mat,b,uv);
+ // check now if the point is inside the triangle
+ if (uv[0] >= -1.e-6 && uv[1] >= -1.e-6 &&
+ 1.-uv[0]-uv[1] >= -1.e-6 ) {
+ SVector3 pp =
+ ct->p[0] * ( 1.-uv[0]-uv[1] ) +
+ ct->p[1] * uv[0] +
+ ct->p[2] * uv[1] ;
+ return GPoint(s[IT].x(),s[IT].y(),s[IT].z(),this,pp);
+ }
+ }
+ }
+ }
+ GPoint pp(0);
+ pp.setNoSuccess();
+ Msg::Debug("ARGG no success intersection circle");
+ // Msg::Info("ARGG no success intersection circle");
+ // printf("Point(1) = {%g,%g,%g};\n",p.x(),p.y(),p.z());
+ // printf("Point(2) = {%g,%g,%g};\n",p.x()+d*n1.x(),p.y()+d*n1.y(),p.z()+d*n1.z());
+ // printf("Point(3) = {%g,%g,%g};\n",p.x()+d*n2.x(),p.y()+d*n2.y(),p.z()+d*n2.z());
+
+ // // printf("Circle(4) = {2,1,3};\n");
+ // printf("{%g,%g,%g};\n",n1.x(),n1.y(),n1.z());
+ // printf("{%g,%g,%g};\n",n2.x(),n2.y(),n2.z());
+ // printf("coucou --> \n");
+ // if (allProblems){
+ // fprintf(allProblems,"VP(%g,%g,%g){%g,%g,%g};\n",p.x(),p.y(),p.z(),R*n2.x(),R*n2.y(),R*n2.z());
+ // }
+ // getchar();
+ return pp;
+}
+
+
+GPoint discreteDiskFace::intersectionWithCircle2(const SVector3 &n1, const SVector3 &n2,
const SVector3 &p, const double &d,
double uv[2]) const
{
+ // n2 is exterior
SVector3 n = crossprod(n1,n2);
n.normalize();
for (int i=-1;i<(int)(discrete_triangles.size()-geoTriangulation->fillingHoles.size());i++){
@@ -1039,47 +1122,46 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
const double c = dot(q-p,q-p) - d*d;
const double delta = b*b-4*a*c;
if (delta >= 0){
- const double ta = (-b + sqrt(delta)) / (2.*a);
- const double tb = (-b - sqrt(delta)) / (2.*a);
- SVector3 s1 = q + m * ta;
- SVector3 s2 = q + m * tb;
- SVector3 s;
- if (dot(s1-p,n1) > 0){
- s = s1;
- }
- else if (dot(s2-p,n1) > 0){
- s = s2;
- }
- else continue;
-
- // we have now to look if the point is inside the triangle
- // s = v1 + u t1 + v t2
- // we know the system has a solution because s is in the plane
- // defined by v1 and v2 we solve then
- // (s - v1) . t1 = u t1^2 + v t2 . t1
- // (s - v1) . t2 = u t1 . t2 + v t2^2
-
- double mat[2][2], b[2],uv[2];
- mat[0][0] = dot(t1,t1);
- mat[1][1] = dot(t2,t2);
- mat[0][1] = mat[1][0] = dot(t1,t2);
- b[0] = dot(s-v0,t1) ;
- b[1] = dot(s-v0,t2) ;
- sys2x2(mat,b,uv);
- // check now if the point is inside the triangle
- if (uv[0] >= -1.e-6 && uv[1] >= -1.e-6 &&
- 1.-uv[0]-uv[1] >= -1.e-6 ) {
- SVector3 pp =
- ct->p[0] * ( 1.-uv[0]-uv[1] ) +
- ct->p[1] * uv[0] +
- ct->p[2] * uv[1] ;
- return GPoint(s.x(),s.y(),s.z(),this,pp);
+ double sign = (dot(n2,m) > 0)? 1.0:-1.0;
+ const double ta = (-b + sign*sqrt(delta)) / (2.*a);
+ const double tb = (-b - sign*sqrt(delta)) / (2.*a);
+ SVector3 s[2] = {q + m * ta, q + m * tb};
+ for (int IT=0;IT<2;IT++){
+ double mat[2][2], b[2],uv[2];
+ mat[0][0] = dot(t1,t1);
+ mat[1][1] = dot(t2,t2);
+ mat[0][1] = mat[1][0] = dot(t1,t2);
+ b[0] = dot(s[IT]-v0,t1) ;
+ b[1] = dot(s[IT]-v0,t2) ;
+ sys2x2(mat,b,uv);
+ // check now if the point is inside the triangle
+ if (uv[0] >= -1.e-6 && uv[1] >= -1.e-6 &&
+ 1.-uv[0]-uv[1] >= -1.e-6 ) {
+ SVector3 pp =
+ ct->p[0] * ( 1.-uv[0]-uv[1] ) +
+ ct->p[1] * uv[0] +
+ ct->p[2] * uv[1] ;
+ return GPoint(s[IT].x(),s[IT].y(),s[IT].z(),this,pp);
+ }
}
}
}
GPoint pp(0);
pp.setNoSuccess();
- Msg::Debug("ARGG no success intersection circle");
+ // Msg::Debug("ARGG no success intersection circle");
+ Msg::Info("ARGG no success intersection circle");
+ // printf("Point(1) = {%g,%g,%g};\n",p.x(),p.y(),p.z());
+ // printf("Point(2) = {%g,%g,%g};\n",p.x()+d*n1.x(),p.y()+d*n1.y(),p.z()+d*n1.z());
+ // printf("Point(3) = {%g,%g,%g};\n",p.x()+d*n2.x(),p.y()+d*n2.y(),p.z()+d*n2.z());
+
+ // // printf("Circle(4) = {2,1,3};\n");
+ // printf("{%g,%g,%g};\n",n1.x(),n1.y(),n1.z());
+ // printf("{%g,%g,%g};\n",n2.x(),n2.y(),n2.z());
+ // printf("coucou --> \n");
+ // if (allProblems){
+ // fprintf(allProblems,"VP(%g,%g,%g){%g,%g,%g};\n",p.x(),p.y(),p.z(),d*n2.x(),d*n2.y(),d*n2.z());
+ // }
+ // getchar();
return pp;
}
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index d0940c62b1f9fa4a36e1dd0a2bda5c6c2ad2015c..163419cfe30bd770d16ff71d7d31efaccdc4bfb0 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -221,9 +221,8 @@ class discreteDiskFaceTriangle {
class discreteDiskFace : public GFace {
GFace *_parent;
void buildOct(std::vector<GFace*> *CAD = NULL) const;
- bool parametrize(bool);// const;
- void checklsys(linearSystemCSRTaucs<double>*,dofManager<double>*,int);
- void putOnView(bool,bool);
+ bool parametrize();// const;
+ void checklsys(linearSystemCSR<double>*,dofManager<double>*,int);
bool checkOrientationUV();
void optimize();
@@ -245,8 +244,12 @@ class discreteDiskFace : public GFace {
GPoint intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
const SVector3 &p, const double &d,
double uv[2]) const;
+ GPoint intersectionWithCircle2(const SVector3 &n1, const SVector3 &n2,
+ const SVector3 &p, const double &d,
+ double uv[2]) const;
void printAtlasMesh () ;
void printParamMesh () ;
+ void putOnView(int iFace, int iMap, bool,bool);
std::vector<MElement*> discrete_triangles;
protected:
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 4e90429bf8a7948725d61c7f38cca4dc7d4611b6..795ac1a6b3e73b110d6926d4f3902dc9a1e2e1c3 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -170,7 +170,8 @@ void discreteFace::createGeometry()
for(unsigned int i=0; i<toParam.size(); i++){
discreteDiskFace *df = new discreteDiskFace (this,toParam[i], order,(_CAD.empty() ? NULL : &_CAD));
- df->printAtlasMesh();
+ // df->printAtlasMesh();
+ // df->putOnView(tag(), i, true,true);
df->replaceEdges(toParam[i]->my_GEdges);
_atlas.push_back(df);
}
@@ -195,11 +196,14 @@ void discreteFace::gatherMeshes()
SPoint2 pc = (p0+p1+p2)*(1./3.0);
discreteDiskFaceTriangle *mt=NULL;
double xi, eta;
+ // FIXME THAT SUCKS !!!!!
_atlas[i]->getTriangleUV(pc.x(),pc.y(), &mt, xi,eta);
if (mt && mt->gf)mt->gf->triangles.push_back(t);
- else Msg::Warning ("FILE %s LINE %d Triangle from atlas part %u has no classification for (u;v)=(%f;%f)",__FILE__,__LINE__,i+1,pc.x(),pc.y());
-
-
+ else {
+ if (_CAD.empty())triangles.push_back(t);
+ else
+ Msg::Warning ("FILE %s LINE %d Triangle from atlas part %u has no classification for (u;v)=(%f;%f)",__FILE__,__LINE__,i+1,pc.x(),pc.y());
+ }
}
_atlas[i]->triangles.clear();
@@ -226,13 +230,22 @@ void discreteFace::mesh(bool verbose)
#if defined(HAVE_ANN) && defined(HAVE_SOLVER) && defined(HAVE_MESH)
if (!CTX::instance()->meshDiscrete) return;
+ Msg::Info("Discrete Face %d is going to be meshed",tag());
for (unsigned int i=0;i<_atlas.size();i++){
- Msg::Info("Discrete Face %d is going to be meshed",i);
- _atlas[i]->mesh(verbose);/*
- const char *name = "atlas%d";
- char filename[256];
- sprintf(filename,name,i);t0
- _atlas[i]->printPhysicalMesh(filename);*/
+ // {
+ // void openProblemsON(void);
+ // if (tag() == 3) openProblemsON();
+ // }
+ _atlas[i]->mesh(verbose);
+ // {
+ // void openProblemsOFF(void);
+ // if (tag()==3) openProblemsOFF();
+ // }
+ /*
+ const char *name = "atlas%d";
+ char filename[256];
+ sprintf(filename,name,i);t0
+ _atlas[i]->printPhysicalMesh(filename);*/
}
gatherMeshes();
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index f2c8ca206edd4afcee68ac2a10403427bebbb45c..c0cd5284b5bd3a291beab688fbb583d242bcf305 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -745,11 +745,9 @@ bool insertVertexB (std::list<edgeXface> &shell,
double d3 = distance(it->v[0],it->v[1]);
// avoid angles that are too obtuse
- double cosv = ((d1*d1+d2*d2-d3*d3)/(2.*d1*d2));
-
- if ((d1 < LL * .25 || d2 < LL * .25 || cosv < -.9999) && !force) {
+ if ((d1 < LL * .5 || d2 < LL * .5 || ((d1*d1+d2*d2-d3*d3)/(2.*d1*d2)) < -.9999) && !force) {
onePointIsTooClose = true;
- // printf("%12.5E %12.5E %12.5E %12.5E \n",d1,d2,LL,cosv);
+ // printf("%12.5E %12.5E %12.5E \n",d1,d2,LL);
}
newTris[k++] = t4;
@@ -1290,38 +1288,39 @@ bool optimalPointFrontalB (GFace *gf,
{
// as a starting point, let us use the "fast algo"
double d = optimalPointFrontal (gf,worst,active_edge,data,newPoint,metric);
- int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
+ int ip1 = (active_edge + 2)%3;
int ip2 = active_edge;
+ int ip3 = (active_edge + 1)%3;
MVertex *v1 = worst->tri()->getVertex(ip1);
MVertex *v2 = worst->tri()->getVertex(ip2);
- MTriangle *t = worst->tri();
- double p1[3] = {t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()};
- double p2[3] = {t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()};
- double p3[3] = {t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z()};
- double c[3];
- circumCenterXYZ(p1, p2, p3, c);
+ MVertex *v3 = worst->tri()->getVertex(ip3);
SVector3 middle ((v1->x()+v2->x())*.5,(v1->y()+v2->y())*.5,(v1->z()+v2->z())*.5);
- SVector3 center(c[0],c[1],c[2]);
- SVector3 v1v2 (v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
- SVector3 n1 = center - middle;
- SVector3 n2 = crossprod(v1v2,n1);
+ SVector3 v1v2 (v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
+ SVector3 tmp (v3->x()-middle.x(),v3->y()-middle.y(),v3->z()-middle.z());
+
+ SVector3 n1 = crossprod(v1v2,tmp);
+ SVector3 n2 = crossprod(n1,v1v2); // n2 is exterior
n1.normalize();
n2.normalize();
- // we look for a point that is
- // P = d * (n1 cos(t) + n2 sin(t)) that is on the surface
- // so we have to find t, starting with t = 0
#if defined(HAVE_ANN) && defined(HAVE_SOLVER)
if (gf->geomType() == GEntity::DiscreteDiskSurface){
discreteDiskFace *ddf = dynamic_cast<discreteDiskFace*> (gf);
if (ddf){
+ // not optimal @ all !!!
GPoint gp = ddf->intersectionWithCircle(n1,n2,middle,d,newPoint);
if (gp.succeeded()){
newPoint[0] = gp.u();
newPoint[1] = gp.v();
return true;
}
+ // gp = ddf->intersectionWithCircle(n2,n1,middle,d,newPoint);
+ // if (gp.succeeded()){
+ // newPoint[0] = gp.u();
+ // newPoint[1] = gp.v();
+ // return true;
+ // }
return false;
}
}
@@ -1359,6 +1358,7 @@ void bowyerWatsonFrontal(GFace *gf,
std::map<MVertex* , MVertex*>* equivalence,
std::map<MVertex*, SPoint2> * parametricCoordinates)
{
+
std::set<MTri3*,compareTri3Ptr> AllTris;
std::set<MTri3*,compareTri3Ptr> ActiveTris;
bidimMeshData DATA(equivalence,parametricCoordinates);
@@ -1382,21 +1382,13 @@ void bowyerWatsonFrontal(GFace *gf,
int ITERATION = 0;
while (1){
++ITERATION;
- /* if(ITERATION % 10== 0 && CTX::instance()->mesh.saveAll){
- char name[245];
- sprintf(name,"delFrontal_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
- _printTris (name, AllTris.begin(), AllTris.end(), &DATA);
- sprintf(name,"delFrontal_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
- _printTris (name, ActiveTris.begin(), ActiveTris.end(), &DATA);
- }*/
- /* if(ITER % 100== 0){
- char name[245];
- sprintf(name,"delfr2d%d-ITER%d.pos",gf->tag(),ITER);
- _printTris (name, AllTris, Us,Vs,true);
- // sprintf(name,"delfr2dA%d-ITER%d.pos",gf->tag(),ITER);
- // _printTris (name, ActiveTris, Us,Vs,false);
- }
- */
+ if(ITERATION % 10== 0 && CTX::instance()->mesh.saveAll){
+ char name[245];
+ sprintf(name,"delFrontal_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
+ _printTris (name, AllTris.begin(), AllTris.end(), NULL);
+ sprintf(name,"delFrontal_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
+ _printTris (name, ActiveTris.begin(), ActiveTris.end(), NULL);
+ }
// printf("active_tris.size = %d\n",ActiveTris.size());
if (!ActiveTris.size())break;
MTri3 *worst = (*ActiveTris.begin());
@@ -1420,8 +1412,6 @@ void bowyerWatsonFrontal(GFace *gf,
}
*/
}
-
-
nbSwaps = edgeSwapPass(gf, AllTris, SWCR_QUAL, DATA);
/*
char name[245];
@@ -1447,6 +1437,8 @@ void bowyerWatsonFrontal(GFace *gf,
}
}
#endif
+ // getchar();
+
}
void optimalPointFrontalQuad (GFace *gf,
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 6b7b00512cd794b31ebcfc406d5ea4671c9148dd..a6b1963883f8859b595cf8d2b04e84f21e417010 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -836,6 +836,11 @@ void getAllBoundaryLayerVertices (GFace *gf, std::set<MVertex*> &vs){
void laplaceSmoothing(GFace *gf, int niter, bool infinity_norm)
{
if (!niter)return;
+ RelocateVertices (gf, niter);
+ return;
+
+
+ printf("argh\n");
std::set<MVertex*> vs;
getAllBoundaryLayerVertices (gf, vs);
v2t_cont adj;
@@ -872,6 +877,12 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
MTri3 *t2 = t1->getNeigh(iLocalEdge);
if(!t2) return false;
+ double triQualityRef = 0;
+ if (cr == SWCR_QUAL){
+ triQualityRef=std::min(qmTriangle::gamma(t1->tri()), qmTriangle::gamma(t2->tri()));
+ if (triQualityRef > 0.5)return false;
+ }
+
MVertex *v1 = t1->tri()->getVertex(iLocalEdge == 0 ? 2 : iLocalEdge - 1);
MVertex *v2 = t1->tri()->getVertex((iLocalEdge) % 3);
MVertex *v3 = t1->tri()->getVertex((iLocalEdge + 1) % 3);
@@ -882,12 +893,14 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
for(int i = 0; i < 3; i++)
- if(t2->tri()->getVertex(i) != v1 && t2->tri()->getVertex(i) != v2)
+ if(t2->tri()->getVertex(i) != v1 && t2->tri()->getVertex(i) != v2){
v4 = t2->tri()->getVertex(i);
+ break;
+ }
- swapquad sq (v1, v2, v3, v4);
- if(configs.find(sq) != configs.end()) return false;
- configs.insert(sq);
+ // swapquad sq (v1, v2, v3, v4);
+ // if(configs.find(sq) != configs.end()) return false;
+ // configs.insert(sq);
if (edgeSwapDelProj(v3,v4,v2,v1))return false;
@@ -895,15 +908,15 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
MTriangle *t1b = new MTriangle(v2, v3, v4);
MTriangle *t2b = new MTriangle(v4, v3, v1);
+
switch(cr){
case SWCR_QUAL:
{
- const double triQualityRef = std::min(qmTriangle::gamma(t1->tri()),
- qmTriangle::gamma(t2->tri()));
const double triQuality = std::min(qmTriangle::gamma(t1b),
qmTriangle::gamma(t2b));
+
if (!edgeSwapDelProj(v1,v2,v3,v4)){
- if(triQuality < triQualityRef){
+ if(triQuality < triQualityRef/* && triQuality < .001*/){
delete t1b;
delete t2b;
return false;
@@ -999,7 +1012,7 @@ int edgeSwapPass(GFace *gf, std::set<MTri3*, compareTri3Ptr> &allTris,
int nbSwapTot = 0;
std::set<swapquad> configs;
- for(int iter = 0; iter < 1200; iter++){
+ for(int iter = 0; iter < 3; iter++){
int nbSwap = 0;
std::vector<MTri3*> newTris;
for(CONTAINER::iterator it = allTris.begin(); it != allTris.end(); ++it){
@@ -1011,17 +1024,12 @@ int edgeSwapPass(GFace *gf, std::set<MTri3*, compareTri3Ptr> &allTris,
}
}
}
- else{
- delete *it;
- CONTAINER::iterator itb = it;
- ++it;
- allTris.erase(itb);
- if(it == allTris.end()) break;
- }
}
+
allTris.insert(newTris.begin(), newTris.end());
nbSwapTot += nbSwap;
if(nbSwap == 0) break;
+ // printf("%d swaps %d new tris\n",nbSwap, newTris.size());
}
return nbSwapTot;
}
diff --git a/Mesh/meshGRegionRelocateVertex.cpp b/Mesh/meshGRegionRelocateVertex.cpp
index 0536ceb23e8d408d101081bc06356f0f92fc5e38..b7b7185eed13d11a8c6df6b89bd27d79ba4ea119 100644
--- a/Mesh/meshGRegionRelocateVertex.cpp
+++ b/Mesh/meshGRegionRelocateVertex.cpp
@@ -212,8 +212,10 @@ static double _relocateVertex(GFace* gf, MVertex *ver,
p1 *= 1./(double)counter;
double worst;
double xi = Maximize_Quality_Golden_Section( ver, gf, p1, p2, lt , tol, worst);
+ // if (xi != 0)printf("xi = %g\n",xi);
SPoint2 p = p1*(1-xi) + p2*xi;
GPoint pp = gf->point(p);
+ if (!pp.succeeded())return 2.0;
ver->x() = pp.x();
ver->y() = pp.y();
ver->z() = pp.z();
diff --git a/Mesh/meshRefine.cpp b/Mesh/meshRefine.cpp
index 11ba91b6c747a53d581f275c08a1d3cebaf34540..e4255ff9b5d76655df36d4d53ee0814e49ed7786 100644
--- a/Mesh/meshRefine.cpp
+++ b/Mesh/meshRefine.cpp
@@ -449,6 +449,9 @@ static void Subdivide(GRegion *gr, bool splitIntoHexas, faceContainer &faceVerti
void RefineMesh(GModel *m, bool linear, bool splitIntoQuads, bool splitIntoHexas)
{
+ // splitIntoQuads = true;
+ // splitIntoHexas = true;
+
Msg::StatusBar(true, "Refining mesh...");
double t1 = Cpu();
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index 37f92fa92da480739be59214933ec39e6337d756..da08e267842682aee73ab6174f4daf3896b9d2bb 100644
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -422,7 +422,7 @@ namespace {
template <class T>
bool find_value_in_multiset(const std::multiset<T>& set, const T& value, T& value_in_set) {
- std::multiset<T>::const_iterator it = set.find(value);
+ typename std::multiset<T>::const_iterator it = set.find(value);
for (; it != set.end() && it->get_hash() == value.get_hash(); ++it) {
if (value.same_vertices(*it)) {
value_in_set = *it;
@@ -435,7 +435,7 @@ namespace {
// Check if one value exist in a multiset hash table of Recombinator
template <class T>
bool find_value_in_multiset(const std::multiset<T>& set, const T& value) {
- std::multiset<T>::const_iterator it = set.find(value);
+ typename std::multiset<T>::const_iterator it = set.find(value);
for (; it != set.end() && it->get_hash() == value.get_hash(); ++it) {
if (value.same_vertices(*it)) {
return true;
@@ -5785,7 +5785,7 @@ bool cliques_losses_graph<T>::compatibility(
const T &u, const hash_key &u_key,
const T &v, const hash_key &v_key)
{
- return !cliques_compatibility_graph::compatibility(u, u_key, v, v_key);
+ return !cliques_compatibility_graph<T>::compatibility(u, u_key, v, v_key);
}
diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index a2dca415acedbba5c28babc8c781a2060e264b4b..94d620b8485979fde6088f2e714914e9560a2679 100644
--- a/Solver/linearSystem.h
+++ b/Solver/linearSystem.h
@@ -41,7 +41,7 @@ class linearSystem : public linearSystemBase {
virtual ~linearSystem (){}
virtual void addToMatrix(int _row, int _col, const scalar &val) = 0;
virtual void getFromMatrix(int _row, int _col, scalar &val) const = 0;
- virtual void addToRightHandSide(int _row, const scalar &val) = 0;
+ virtual void addToRightHandSide(int _row, const scalar &val, int ith = 0) = 0;
virtual void getFromRightHandSide(int _row, scalar &val) const = 0;
virtual void getFromSolution(int _row, scalar &val) const = 0;
virtual void addToSolution(int _row, const scalar &val) = 0;
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 1e73598062aa88836971e1811a8c8c5d85d9d163..d25cd8dd1e1652942bb2980a000f740f80cdbec8 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/Solver/linearSystemCSR.cpp
@@ -490,6 +490,9 @@ int linearSystemCSRTaucs<double>::systemSolve()
}
sorted = true;
+ // taucs_logfile("stderr");
+
+
taucs_ccs_matrix myVeryCuteTaucsMatrix;
myVeryCuteTaucsMatrix.n = myVeryCuteTaucsMatrix.m = _b->size();
//myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)_ptr->array;
@@ -497,9 +500,13 @@ int linearSystemCSRTaucs<double>::systemSolve()
myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)_ai->array;
myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)_jptr->array;
myVeryCuteTaucsMatrix.values.d = (double*)_a->array;
- myVeryCuteTaucsMatrix.flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
+ if (_symmetric)
+ myVeryCuteTaucsMatrix.flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
+ else
+ myVeryCuteTaucsMatrix.flags = TAUCS_DOUBLE ;
- char* options[] = { (char*)"taucs.factor.LLT=true", NULL };
+ char* options[] = { _symmetric ? (char*)"taucs.factor.LLT=true" : (char*)"taucs.factor.LLT=true",NULL };
+ // char* options[] = { NULL };
double t1 = Cpu();
int result = taucs_linsolve(&myVeryCuteTaucsMatrix,
NULL,
diff --git a/Solver/linearSystemCSR.h b/Solver/linearSystemCSR.h
index 17eb46e9a5630d9e02c87ae75aadd0abfbf56292..93e056dd62e70f45e549e5e32a8b441619f0b8ee 100644
--- a/Solver/linearSystemCSR.h
+++ b/Solver/linearSystemCSR.h
@@ -35,6 +35,8 @@ class linearSystemCSR : public linearSystem<scalar> {
std::vector<scalar> *_b, *_x;
sparsityPattern _sparsity; // only used for pre-allocation, does not store the sparsity once allocated
public:
+ int getNNZ() {return CSRList_Nbr(linearSystemCSR<scalar>::_a);}
+ int getNbUnk() {return linearSystemCSR<scalar>::_b->size();}
linearSystemCSR()
: sorted(false), _entriesPreAllocated(false), _a(0), _b(0), _x(0) {}
virtual bool isAllocated() const { return _a != 0; }
@@ -121,7 +123,7 @@ class linearSystemCSR : public linearSystem<scalar> {
{
Msg::Error("getFromMatrix not implemented for CSR");
}
- virtual void addToRightHandSide(int row, const scalar &val)
+ virtual void addToRightHandSide(int row, const scalar &val, int ith = 0)
{
if (!_b) return;
if(val != scalar()) (*_b)[row] += val;
@@ -194,17 +196,18 @@ class linearSystemCSRGmm : public linearSystemCSR<scalar> {
template <class scalar>
class linearSystemCSRTaucs : public linearSystemCSR<scalar> {
+ bool _symmetric;
public:
- linearSystemCSRTaucs(){}
+ linearSystemCSRTaucs(bool s = true) : _symmetric(s){}
virtual ~linearSystemCSRTaucs(){}
virtual void addToMatrix(int il, int ic, const scalar &val)
{
- if (il <= ic) {
+ if (!_symmetric || il <= ic) {
linearSystemCSR<scalar>::addToMatrix(il, ic, val);
}
}
virtual void insertInSparsityPattern(int il, int ic) {
- if (il <= ic)
+ if (!_symmetric || il <= ic)
linearSystemCSR<scalar>::insertInSparsityPattern(il,ic);
}
virtual int systemSolve()
@@ -215,8 +218,6 @@ class linearSystemCSRTaucs : public linearSystemCSR<scalar> {
}
#endif
;
- int getNNZ() {return CSRList_Nbr(linearSystemCSR<scalar>::_a);}
- int getNbUnk() {return linearSystemCSR<scalar>::_b->size();}
};
#endif
diff --git a/Solver/linearSystemFull.h b/Solver/linearSystemFull.h
index 173245dfa5cfdd502e0146bb8d94e08ca1702a8b..1ae37db5fa43495d3ec58a3b748bd983e8616ad1 100644
--- a/Solver/linearSystemFull.h
+++ b/Solver/linearSystemFull.h
@@ -50,7 +50,7 @@ class linearSystemFull : public linearSystem<scalar> {
{
val = (*_a)(row, col);
}
- virtual void addToRightHandSide(int row, const scalar &val)
+ virtual void addToRightHandSide(int row, const scalar &val, int ith=0)
{
if(val != 0.0) (*_b)(row) += val;
}
diff --git a/Solver/linearSystemGMM.h b/Solver/linearSystemGMM.h
index 41da6e76aefcd8e8fe44f0d710124049147418c2..b69092e46686dcb01e42a2cdcdf2566a070e1035 100644
--- a/Solver/linearSystemGMM.h
+++ b/Solver/linearSystemGMM.h
@@ -58,7 +58,7 @@ class linearSystemGmm : public linearSystem<scalar> {
val = (*_a)(row, col);
}
- virtual void addToRightHandSide(int row, const scalar &val)
+ virtual void addToRightHandSide(int row, const scalar &val, int ith = 0)
{
if(val != 0.0) (*_b)[row] += val;
}
@@ -133,7 +133,7 @@ class linearSystemGmm : public linearSystem<scalar> {
virtual void allocate(int nbRows) {}
virtual void addToMatrix(int row, int col, const scalar &val) {}
virtual void getFromMatrix(int row, int col, scalar &val) const {}
- virtual void addToRightHandSide(int row, const scalar &val) {}
+ virtual void addToRightHandSide(int row, const scalar &val, int ith = 0) {}
virtual void getFromRightHandSide(int row, scalar &val) const {}
virtual void addToSolution(int row, const scalar &val) {}
virtual void getFromSolution(int row, scalar &val) const {}
diff --git a/Solver/linearSystemMUMPS.cpp b/Solver/linearSystemMUMPS.cpp
index fe09aae114e139da9f7650a6a0d5b9969fc80792..1f6349d5b74d8fc51c2a9fb4790cf0a4ef377d1a 100644
--- a/Solver/linearSystemMUMPS.cpp
+++ b/Solver/linearSystemMUMPS.cpp
@@ -53,8 +53,8 @@ bool linearSystemMUMPS<double>::isAllocated() const
void linearSystemMUMPS<double>::allocate(int nbRows)
{
_n = nbRows;
- _b.reserve(_n);
- _x.reserve(_n);
+ _b.resize(_n);
+ _x.resize(_n);
_a.reserve(10*_n);
_irn.reserve(10*_n);
_jcn.reserve(10*_n);
@@ -107,7 +107,7 @@ int linearSystemMUMPS<double>::systemSolve()
id.comm_fortran = USE_COMM_WORLD;
- Msg::Info("MUMPS initialization");
+ Msg::Debug("MUMPS initialization");
id.job=-1;
dmumps_c(&id);
mumpserror(id.info[0], id.info[1]);
@@ -130,15 +130,15 @@ int linearSystemMUMPS<double>::systemSolve()
id.icntl[5-1]=0;
id.icntl[18-1]=0;
- Msg::Info("MUMPS analysis, LU factorization, and back substitution");
+ Msg::Debug("MUMPS analysis, LU factorization, and back substitution");
id.job = 6;
dmumps_c(&id);
mumpserror(id.info[0], id.info[1]);
- Msg::Info("MUMPS destroy");
+ Msg::Debug("MUMPS destroy");
id.job=-2;
dmumps_c(&id);
- Msg::Info("MUMPS end");
+ Msg::Debug("MUMPS end");
mumpserror(id.info[0], id.info[1]);
_x.clear();
@@ -216,8 +216,9 @@ void linearSystemMUMPS<double>::getFromMatrix(int row, int col, double &val) con
else val = _a[it->second];
}
-void linearSystemMUMPS<double>::addToRightHandSide(int row, const double &val)
+void linearSystemMUMPS<double>::addToRightHandSide(int row, const double &val, int ith)
{
+ // printf("adding %g to %d\n",val,row);
if((int)_b.size() <= row) {
_b.resize(row+1);
_b[row] = val;
@@ -230,11 +231,12 @@ void linearSystemMUMPS<double>::addToRightHandSide(int row, const double &val)
void linearSystemMUMPS<double>::getFromRightHandSide(int row, double &val) const
{
if((int)_b.size() <= row) val = 0.;
- else val = _b[row];
+ val = _b[row];
}
void linearSystemMUMPS<double>::getFromSolution(int row, double &val) const
{
+ // printf("x[%d] = %g\n",row,_x[row]);
if((int)_x.size() <= row) val = 0.;
else val = _x[row];
}
diff --git a/Solver/linearSystemMUMPS.h b/Solver/linearSystemMUMPS.h
index e8b953a5bb7ee8475450e3ffd9925d767b85bd25..0b6cd87f2293b56e888ed3ad0731f2d0b06d37be 100644
--- a/Solver/linearSystemMUMPS.h
+++ b/Solver/linearSystemMUMPS.h
@@ -38,7 +38,7 @@ class linearSystemMUMPS : public linearSystem<scalar> {
virtual void addToMatrix(int row, int col, const double &val) {}
virtual void getFromMatrix(int row, int col, double &val) const {}
- virtual void addToRightHandSide(int row, const scalar &val) {}
+ virtual void addToRightHandSide(int row, const scalar &val, int ith = 0) {}
virtual void getFromRightHandSide(int row, scalar &val) const {}
virtual void getFromSolution(int row, scalar &val) const {}
virtual void addToSolution(int row, const scalar &val) {}
@@ -80,7 +80,7 @@ class linearSystemMUMPS<double> : public linearSystem<double> {
virtual void addToMatrix(int row, int col, const double &val);
virtual void getFromMatrix(int row, int col, double &val) const;
- virtual void addToRightHandSide(int row, const double &val);
+ virtual void addToRightHandSide(int row, const double &val, int ith = 0);
virtual void getFromRightHandSide(int row, double &val) const;
virtual void getFromSolution(int row, double &val) const;
virtual void addToSolution(int row, const double &val);
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 76db880392ea9190dbac544e149dcd52bbc924ca..7f4bb8b34567bc657f300562f4398c1fb10a46bc 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -77,7 +77,7 @@ void linearSystemPETSc<fullMatrix<double> >::addToMatrix(int row, int col,
template<>
void linearSystemPETSc<fullMatrix<double> >::addToRightHandSide(int row,
- const fullMatrix<double> &val)
+ const fullMatrix<double> &val, int ith)
{
if (!_entriesPreAllocated)
preAllocateEntries();
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index 646e065d075544df325058d6cf4eae3b8827a902..19cc01ad7e0c172701460f2e57eabc09225822ac 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -102,7 +102,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
void print();
void clear();
void getFromMatrix(int row, int col, scalar &val) const;
- void addToRightHandSide(int row, const scalar &val);
+ void addToRightHandSide(int row, const scalar &val, int ith=0);
void getFromRightHandSide(int row, scalar &val) const;
double normInfRightHandSide() const;
int getNumKspIteration() const;
@@ -139,7 +139,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
void clear(){}
void addToMatrix(int row, int col, const scalar &val) {}
void getFromMatrix(int row, int col, scalar &val) const {}
- void addToRightHandSide(int row, const scalar &val) {}
+ void addToRightHandSide(int row, const scalar &val, int ith = 0) {}
void addToSolution(int row, const scalar &val) {}
void getFromRightHandSide(int row, scalar &val) const {}
void getFromSolution(int row, scalar &val) const {}
diff --git a/benchmarks/2d/BL3.geo b/benchmarks/2d/BL3.geo
index d145eceb2df52586f2673f09900d9a904accb0c1..ac8701a126fb089fd2eb733ac12b7d8e690f4b5a 100644
--- a/benchmarks/2d/BL3.geo
+++ b/benchmarks/2d/BL3.geo
@@ -60,7 +60,7 @@ Field[1].hwall_n = 0.001;
//+
Field[1].ratio = 1.4;
//+
-Field[1].thickness = .15;
+Field[1].thickness = .35;
BoundaryLayer Field = 1;
//+
Line Loop(7) = {1, 2, 3, 4};