diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 5dab88d033c11eab362d792b233511fe63765ec2..2591aa0a3bd0e5276a850033b0a1080646f8019f 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -368,17 +368,19 @@ bool GFaceCompound::checkFolding(std::vector<MVertex*> &ordered) const
int maxSize = (i==0) ? ordered.size()-2: ordered.size()-1;
for(int k = i+2; k < maxSize; ++k){
SPoint3 q1 = coordinates[ordered[k]];
- SPoint3 q2 = coordinates[ordered[k]];
+ SPoint3 q2 = coordinates[ordered[k+1]];
double x[2];
int inters = intersection_segments (p1,p2,q1,q2,x);
if (inters > 0) has_no_folding = false;
}
}
- if ( !has_no_folding )
+ if ( !has_no_folding ) {
Msg::Warning("$$$ Folding for compound face %d", this->tag());
+ }
return has_no_folding;
+
}
//check if the discrete harmonic map is correct
@@ -800,7 +802,6 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
_lsys = _lsysb;
#elif defined(HAVE_TAUCS)
_lsys = new linearSystemCSRTaucs<double>;
- printf("compound with taucs\n");
#else
_lsys = new linearSystemFull<double>;
#endif
@@ -1115,7 +1116,6 @@ bool GFaceCompound::parametrize_conformal_spectral() const
{
#if defined(HAVE_PETSC)
-
std::vector<MVertex*> ordered;
std::vector<double> coords;
bool success = orderVertices(_U0, ordered, coords);
@@ -1159,12 +1159,27 @@ bool GFaceCompound::parametrize_conformal_spectral() const
//-------------------------------
myAssembler.setCurrentMatrix("B");
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
- MVertex *v = *itv;
- myAssembler.numberVertex(v, 0, 1);
- myAssembler.numberVertex(v, 0, 2);
- }
+// massTerm mass1(model(), 1, &MONE);
+// massTerm mass2(model(), 2, &MONE);
+// it = _compound.begin();
+// for( ; it != _compound.end() ; ++it){
+// for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+// SElement se((*it)->triangles[i]);
+// mass1.addToMatrix(myAssembler, &se);
+// mass2.addToMatrix(myAssembler, &se);
+// }
+// }
+
+// std::list<GEdge*>::const_iterator it0 = _U0.begin();
+// for( ; it0 != _U0.end(); ++it0 ){
+// for(unsigned int i = 0; i < (*it0)->lines.size(); i++ ){
+// SElement se((*it0)->lines[i]);
+// mass1.addToMatrix(myAssembler, &se);
+// mass2.addToMatrix(myAssembler, &se);
+// }
+// }
+
double small = 0.0;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
@@ -1176,60 +1191,58 @@ bool GFaceCompound::parametrize_conformal_spectral() const
myAssembler.assemble(v, 0, 1, v, 0, 1, 1.0);
myAssembler.assemble(v, 0, 2, v, 0, 2, 1.0);
}
- }
+ }
-// int NB = ordered.size();
-// for (int i = 0; i< NB; i++){
-// MVertex *v1 = ordered[i];
-// for (int j = i; j< NB; j++){
-// MVertex *v2 = ordered[j];
-// myAssembler.assemble(v1, 0, 1, v2, 0, 1, -1./NB);
-// myAssembler.assemble(v1, 0, 2, v2, 0, 2, -1./NB);
-// myAssembler.assemble(v2, 0, 1, v1, 0, 1, -1./NB);
-// myAssembler.assemble(v2, 0, 2, v1, 0, 2, -1./NB);
-// }
-// }
-
-// diagBCTerm diag1(0, 1, &ONE);
-// diagBCTerm diag2(0, 2, &ONE);
-// it = _compound.begin();
-// for( ; it != _compound.end() ; ++it){
-// for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
-// SElement se((*it)->triangles[i]);
-// diag1.addToMatrix(myAssembler, &se);
-// diag2.addToMatrix(myAssembler, &se);
-// }
-// }
+ int NB = ordered.size();
+ for (int i = 0; i< NB; i++){
+ MVertex *v1 = ordered[i];
+ for (int j = 0; j< NB; j++){
+ MVertex *v2 = ordered[j];
+ myAssembler.assemble(v1, 0, 1, v2, 0, 1, -1./NB);
+ myAssembler.assemble(v1, 0, 2, v2, 0, 2, -1./NB);
+ }
+ }
//-------------------------------
eigenSolver eig(&myAssembler, "B" , "A", true);
- bool converged = eig.solve(1, "largest");
-
+ int nb = 1;
+ bool converged = eig.solve(nb, "largest");
+
if(converged) {
+ double Linfty = 0.0;
int k = 0;
std::vector<std::complex<double> > &ev = eig.getEigenVector(0);
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
- MVertex *v = *itv;
double paramu = ev[k].real();
double paramv = ev[k+1].real();
+ Linfty=std::max(Linfty, sqrt(paramu*paramu+paramv*paramv));
+ k = k+2;
+ }
+ k = 0;
+ double norm2 = 0.0;
+ for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ MVertex *v = *itv;
+ double paramu = ev[k].real()/Linfty;
+ double paramv = ev[k+1].real()/Linfty;
coordinates[v] = SPoint3(paramu,paramv,0.0);
k = k+2;
}
- //if folding take second sallest eigenvalue
-// bool noFolding = checkFolding(ordered);
-// if (!noFolding ){
-// coordinates.clear();
-// int k = 0;
-// std::vector<std::complex<double> > &ev = eig.getEigenVector(1);
-// for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
-// MVertex *v = *itv;
-// double paramu = ev[k].real();
-// double paramv = ev[k+1].real();
-// coordinates[v] = SPoint3(paramu,paramv,0.0);
-// k = k+2;
-// }
-// }
+
+ //if folding take second smallest eigenvalue
+ bool noFolding = checkFolding(ordered);
+ if (!noFolding && nb > 1){
+ coordinates.clear();
+ int k = 0;
+ std::vector<std::complex<double> > &ev = eig.getEigenVector(nb-1);
+ for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ MVertex *v = *itv;
+ double paramu = ev[k].real();
+ double paramv = ev[k+1].real();
+ coordinates[v] = SPoint3(paramu,paramv,0.0);
+ k = k+2;
+ }
+ }
lsysA->clear();
lsysB->clear();
@@ -1238,7 +1251,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
}
else return false;
-
+
#else
return false;
@@ -2046,7 +2059,7 @@ void GFaceCompound::printStuff() const
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
it0->second.y(),it1->second.y(),it2->second.y());
- fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(xyzu,"ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E){%22.15E,%22.15E,%22.15E};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
@@ -2074,7 +2087,7 @@ void GFaceCompound::printStuff() const
it1->second.x(), it1->second.y(), 0.0,
it2->second.x(), it2->second.y(), 0.0,
area, area, area);
- fprintf(uvx,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(uvx,"ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E){%22.15E,%22.15E,%22.15E};\n",
it0->second.x(), it0->second.y(), 0.0,
it1->second.x(), it1->second.y(), 0.0,
it2->second.x(), it2->second.y(), 0.0,
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index a260641cae21cd0a50f7fe336b02b5cdbd0dcfe0..353665900f35516fa90abd3a519271a4df29dd13 100644
--- a/Solver/eigenSolver.cpp
+++ b/Solver/eigenSolver.cpp
@@ -25,6 +25,7 @@ eigenSolver::eigenSolver(dofManager<double> *manager, std::string A,
bool eigenSolver::solve(int numEigenValues, std::string which)
{
+
if(!_A) return false;
Mat A = _A->getMatrix();
Mat B = _B ? _B->getMatrix() : PETSC_NULL;
@@ -53,7 +54,8 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
// set some default options
_try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
_try(EPSSetTolerances(eps,1.e-7,50));
- _try(EPSSetType(eps,EPSARNOLDI)); //EPSKRYLOVSCHUR is default
+ //_try(EPSSetType(eps, EPSKRYLOVSCHUR)); //default
+ _try(EPSSetType(eps,EPSARNOLDI));
//_try(EPSSetType(eps,EPSARPACK));
//_try(EPSSetType(eps,EPSPOWER));
@@ -146,7 +148,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
ev[i] = std::complex<double>(tmpr[i], tmpi[i]);
#endif
}
- _eigenVectors.push_back(ev);
+ _eigenVectors.push_back(ev);
}
_try(VecDestroy(xr));
_try(VecDestroy(xi));
diff --git a/Solver/function.h b/Solver/function.h
index 54229be318bae278b4ed2e200ae5cd1ca0b3d18d..77e32267240e7f5653a670528d2ccb93dcb36ab7 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -19,7 +19,7 @@ class binding;
// b2) The algo insert in the dataCacheMap the dataCache it will provide (ie the quadrature points, the solution value,...)
// b3) Then the algo give this dataCacheMap to the law and ask the source term of the law. This source term will be given as a dataCache, i.e. a cached value, a function to update it and a node in the dependency tree. In this example, to create this datacache two new dataCache will be required : one for the "wind" and one for "xyz". So, the source dataCache will ask the dataCacheMap for the dataCache called "wind". If this one does not already exist, the dataCacheMap will ask the function called "wind" to create it. In turn, the constructor of the dataCache associated with the "wind" will ask for the for the dataCache "xyz" wich will be created by the function "xyz" if it does not already exist and so on. The dataCacheMap contain also a special dataCache associated with the current element. During it's construction it's dataCache keep references to all the dataCache needed to compute it. When they request for the dataCache they need, the dataCacheMap is informed and the dependency tree is built.
//
-// c) When iterating over elements, for each element : the Algo will change the dataCache element so that it point to the current element. This will invalidate everything that depend on the current on the current element (in this case, xyz, wind (because it depends on xyz) and the source dataCache (because it depends on wind). So, in this case everything (exept the current element dataCache) is marked as invalid. Now the algo want to access the value of the 'source' dataCache this one is invalid, so it is recomputed. To recompute it, the value of the wind is accessed and it's value is recomputed, and so on. Now, if some other function access the wind or the position (xyz), they are already marked as valid and are not re-computed.
+// c) When iterating over elements, for each element : the Algo will change the dataCache element so that it points to the current element. This will invalidate everything that depend on the current on the current element (in this case, xyz, wind (because it depends on xyz) and the source dataCache (because it depends on wind). So, in this case everything (exept the current element dataCache) is marked as invalid. Now the algo want to access the value of the 'source' dataCache this one is invalid, so it is recomputed. To recompute it, the value of the wind is accessed and it's value is recomputed, and so on. Now, if some other function access the wind or the position (xyz), they are already marked as valid and are not re-computed.
//
// d) after the loop over the elements the dataCacheMap is deleted and it delete all its dataCaches.
//
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index cd2516c68ed3ca3b8351de694b691bdb28af1d21..789e61705858f165dd0100c075320a4f647e9fed 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -92,8 +92,7 @@ class helmholtzTerm : public femTerm<scalar> {
for (int k = 0; k <= j; k++){
m(j, k) += (K * (Grads[j][0] * Grads[k][0] +
Grads[j][1] * Grads[k][1] +
- Grads[j][2] * Grads[k][2]) + A * sf[j] * sf[k]) *
- weightDetJ;
+ Grads[j][2] * Grads[k][2]) + A * sf[j] * sf[k]) * weightDetJ;
}
}
}
diff --git a/Solver/laplaceTerm.h b/Solver/laplaceTerm.h
index a118637e472c4d4b70c0ea7a5ca806bf976b4afa..7122c69334e52483687e78f28f81e0d4e857c27d 100644
--- a/Solver/laplaceTerm.h
+++ b/Solver/laplaceTerm.h
@@ -15,4 +15,12 @@ class laplaceTerm : public helmholtzTerm<double> {
: helmholtzTerm<double>(gm, iField, iField, k, 0) {}
};
+// a \nabla U
+class massTerm : public helmholtzTerm<double> {
+ public:
+ massTerm(GModel *gm, int iField, simpleFunction<double> *a)
+ : helmholtzTerm<double>(gm, iField, iField, 0, a) {}
+};
+
+
#endif
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 4bb7d9b9d6f7c21fc5811bca618b22a464c8fdd4..73637c6874988773357f77d50ec26b078ab3591a 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -826,7 +826,6 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
#if defined(HAVE_TAUCS)
_lsys = new linearSystemCSRTaucs<double>;
- printf("taucs again\n");
#elif defined(HAVE_GMM)
linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
_lsysb->setGmres(1);