conformal_spectral map with slepc

Geo/GFaceCompound.cpp

@@ -25,6 +25,7 @@
 #include "distanceTerm.h"
 #include "crossConfTerm.h"
 #include "convexCombinationTerm.h"
+#include "diagBCTerm.h" 
 #include "linearSystemGMM.h"
 #include "linearSystemCSR.h"
 #include "linearSystemFull.h"
@@ -502,13 +503,14 @@ bool GFaceCompound::parametrize() const
   else if (_mapping == CONFORMAL){
     Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
     fillNeumannBCS();
-    bool withoutFolding = parametrize_conformal() ;
+    bool withoutFolding = parametrize_conformal_spectral() ;
     printStuff();
     if ( withoutFolding == false ){
       Msg::Warning("$$$ Parametrization switched to harmonic map");
       parametrize(ITERU,HARMONIC); 
       parametrize(ITERV,HARMONIC);
     }
+    //exit(1);
   }
   //Distance function
   //-----------------
@@ -784,14 +786,14 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
 {
 
   if (!_lsys) {
-#if defined(HAVE_PETSC)
+#if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
     _lsys = new linearSystemPETSc<double>;
-#elif defined(_HAVE_TAUCS) 
-    _lsys = new linearSystemCSRTaucs<double>;
-#elif defined(HAVE_GMM)
+#elif defined(HAVE_GMM) && !defined(HAVE_TAUCS)
     linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
     _lsysb->setGmres(1);
     _lsys = _lsysb;
+#elif defined(_HAVE_TAUCS) 
+    _lsys = new linearSystemCSRTaucs<double>;
 #else
     _lsys = new linearSystemFull<double>;
 #endif
@@ -924,7 +926,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
     v->getParameter(0,tGlob);
  
     //find compound Edge
-      GEdgeCompound *gec = dynamic_cast<GEdgeCompound*>(v->onWhat());
+    GEdgeCompound *gec = dynamic_cast<GEdgeCompound*>(v->onWhat());
      
     if(gec){
 
@@ -1102,6 +1104,88 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
 
 }
 
+bool GFaceCompound::parametrize_conformal_spectral() const
+{
+
+#if defined(HAVE_PETSC)
+
+  std::vector<MVertex*> ordered;
+  std::vector<double> coords;  
+  bool success = orderVertices(_U0, ordered, coords);
+
+  linearSystem <double> *lsysA  = new linearSystemPETSc<double>;
+  linearSystem <double> *lsysB  = new linearSystemPETSc<double>;
+  dofManager<double> myAssembler(lsysA, lsysB);
+
+  //-------------------------------
+  myAssembler.setCurrentMatrix("A");
+
+  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);
+  }
+
+  simpleFunction<double> ONE(1.0);
+  simpleFunction<double> MONE(-1.0 );
+  laplaceTerm laplace1(model(), 1, &ONE);
+  laplaceTerm laplace2(model(), 2, &ONE);
+  crossConfTerm cross12(model(), 1, 2, &ONE);
+  crossConfTerm cross21(model(), 2, 1, &MONE);
+  std::list<GFace*>::const_iterator it = _compound.begin(); 
+  for( ; it != _compound.end() ; ++it){
+    for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+      SElement se((*it)->triangles[i]);
+      laplace1.addToMatrix(myAssembler, &se);
+      laplace2.addToMatrix(myAssembler, &se);
+      cross12.addToMatrix(myAssembler, &se);
+      cross21.addToMatrix(myAssembler, &se);
+    }
+  }
+
+  //-------------------------------
+   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);
+   }
+   diagBCTerm diag(0, 1, &ONE);
+   it = _compound.begin(); 
+   for( ; it != _compound.end() ; ++it){
+     for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+       SElement se((*it)->triangles[i]);
+       diag.addToMatrix(myAssembler, &se);
+     }
+   }
+   myAssembler.setCurrentMatrix("A");
+  //-------------------------------
+   eigenSolver eig(&myAssembler, "A" ); //, "B");
+   eig.solve(2, "smallest");
+   //printf("num eigenvalues =%d \n", eig.getNumEigenValues());
+   
+   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();
+     coordinates[v] = SPoint3(paramu,paramv,0.0);
+     k = k+2;
+   }
+  
+   lsysA->clear();
+   lsysB->clear();
+   
+   //check for folding
+   return checkFolding(ordered);
+
+#else
+   return false;
+
+#endif
+}
 bool GFaceCompound::parametrize_conformal() const
 {
 
@@ -1115,7 +1199,7 @@ bool GFaceCompound::parametrize_conformal() const
     return false;
   }
 
-   MVertex *v1 = ordered[0];
+   MVertex *v1  = ordered[0];
    MVertex *v2  = ordered[(int)ceil((double)ordered.size()/2.)];
 
 //   MVertex *v2 ;  
@@ -1271,9 +1355,6 @@ void GFaceCompound::computeNormals() const
       MTriangle *t = (*it)->triangles[i];
       t->getJacobian(0, 0, 0, J);
       SVector3 n (J[2][0],J[2][1],J[2][2]);
-      //SVector3 d1(J[0][0], J[0][1], J[0][2]);
-      //SVector3 d2(J[1][0], J[1][1], J[1][2]);
-      //SVector3 n = crossprod(d1, d2);
       n.normalize();
       for(int j = 0; j < 3; j++){
 	std::map<MVertex*, SVector3>::iterator itn = _normals.find(t->getVertex(j));

Geo/GFaceCompound.h

@@ -74,6 +74,7 @@ class GFaceCompound : public GFace {
   void buildAllNodes() const; 
   void parametrize(iterationStep, typeOfMapping) const;
   bool parametrize_conformal() const;
+  bool parametrize_conformal_spectral() const;
   void compute_distance() const;
   bool checkOrientation(int iter) const;
   bool checkFolding(std::vector<MVertex*> &ordered) const;

Mesh/meshGFace.cpp

@@ -36,6 +36,7 @@
 #include "CreateFile.h"
 #include "Context.h"
 #include "multiscalePartition.h"
+#include "meshGFaceLloyd.h"
 
 void fourthPoint(double *p1, double *p2, double *p3, double *p4)
 {
@@ -1403,11 +1404,13 @@ void partitionAndRemesh(GFaceCompound *gf)
 
   Msg::Info("*** Starting Mesh of surface %d ...", gf->tag());
 
+  //lloydAlgorithm lloyd(10,0);
   for (int i=0; i < NF; i++){
     GFace *gfc =  gf->model()->getFaceByTag(numf + NF + i );
     meshGFace mgf;
     mgf(gfc);
-      
+    //lloyd(gfc);
+
     for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
       MTriangle *t = gfc->triangles[j];
       std::vector<MVertex *> v(3);

Mesh/meshGFaceLloyd.cpp

@@ -38,7 +38,11 @@ void lloydAlgorithm::operator () ( GFace * gf) {
   for (std::set<MVertex*>::iterator it = all.begin(); it != all.end(); ++it){
     SPoint2 p;
     bool success = reparamMeshVertexOnFace(*it, gf, p);
-    if (!success) return;
+    if (!success) {
+      printf("loyd no succes exit \n");
+      exit(1);
+      return;
+    }
     double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / 
       (double)RAND_MAX;
     double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / 

Solver/convexCombinationTerm.h

@@ -51,6 +51,21 @@ class convexCombinationTerm : public femTerm<double> {
     }
   }
 
+//diag matrix
+/*   virtual void elementMatrix(SElement *se, fullMatrix<double> &m) const */
+/*   { */
+
+/*     MElement *e = se->getMeshElement(); */
+/*     m.setAll(0.); */
+/*     const int nbNodes = e->getNumVertices(); */
+/*     for (int j = 0; j < nbNodes; j++){ */
+/*       for (int k = 0; k < nbNodes; k++)  m(j,k) = 0.0; */
+/*       MVertex *v = e->getVertex(j); */
+/*       if( v->onWhat()->dim() < 2) m(j,j) = (nbNodes - 1) ; */
+/*       else m(j,j) = 0.0; */
+/*     } */
+/*   } */
+
 
 };
 

Solver/diagBCTerm.h

+
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _DIAGBC_TERM_H_
+#define _DIAGBC_TERM_H_
+
+#include <assert.h>
+#include "femTerm.h"
+#include "simpleFunction.h"
+#include "Gmsh.h"
+#include "GModel.h"
+#include "SElement.h"
+#include "fullMatrix.h"
+#include "Numeric.h"
+
+class diagBCTerm : public femTerm<double> {
+ protected:
+  const simpleFunction<double> *_k;
+  const int _iField;
+ public:
+  diagBCTerm(GModel *gm, int iField, simpleFunction<double> *k)
+    : femTerm<double>(gm), _iField(iField), _k(k) {}
+  virtual int sizeOfR(SElement *se) const
+  {
+    return se->getMeshElement()->getNumVertices();
+  }
+  virtual int sizeOfC(SElement *se) const 
+  { 
+    return se->getMeshElement()->getNumVertices();
+  }
+  Dof getLocalDofR(SElement *se, int iRow) const
+  {
+    return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), 
+               Dof::createTypeWithTwoInts(0, _iField));
+  }
+
+  virtual void elementMatrix(SElement *se, fullMatrix<double> &m) const
+  {
+
+    MElement *e = se->getMeshElement();
+    m.setAll(0.);
+    const int nbNodes = e->getNumVertices();
+    for (int j = 0; j < nbNodes; j++){
+      for (int k = 0; k < nbNodes; k++) {
+	m(j,k) = 0.0;
+	m(k,j) = 0.0;
+      }
+      m(j,j)  = 1.0;
+/*       MVertex *v = e->getVertex(j); */
+/*       if( v->onWhat()->dim() < 2 )  m(j,j) = (nbNodes - 1) ; */
+/*       else m(j,j) = 0.0; */
+    }
+  } 
+  
+
+};
+
+#endif

Solver/dofManager.h

@@ -101,7 +101,13 @@ class dofManager{
   linearSystem<dataMat> *_current;
    
  public:
-  dofManager(linearSystem<dataMat> *l) : _current(l) { _linearSystems["A"] = l; }
+ dofManager(linearSystem<dataMat> *l) : _current(l) { _linearSystems["A"] = l; }
+ dofManager(linearSystem<dataMat> *l1, linearSystem<dataMat> *l2) : _current(l1) { 
+    _linearSystems.insert(std::make_pair("A", l1));
+    _linearSystems.insert(std::make_pair("B", l2));
+    //_linearSystems.["A"] = l1; 
+    //_linearSystems["B"] = l2;
+  }
   inline void fixDof(Dof key, const dataVec &value)
   {
     fixed[key] = value;
@@ -312,6 +318,15 @@ class dofManager{
     _current->zeroMatrix();
     _current->zeroRightHandSide();
   }  
+  inline void setCurrentMatrix(std::string name){
+    typename std::map<const std::string, linearSystem<dataMat>*>::iterator it =  _linearSystems.find(name);
+     if(it != _linearSystems.end())
+       _current = it->second;
+     else{
+       Msg::Error("Current matrix %s not found ",  name.c_str());
+       throw;
+     }
+  }
   linearSystem<dataMat> *getLinearSystem(std::string &name)
   {
     typename std::map<const std::string, linearSystem<dataMat>*>::iterator it =

Solver/eigenSolver.cpp

@@ -27,12 +27,12 @@ void eigenSolver::solve(int numEigenValues, std::string which)
   if(!_A) return;
   Mat A = _A->getMatrix();
   Mat B = _B ? _B->getMatrix() : PETSC_NULL;
-
+   
   _try(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
   _try(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
   PetscInt N, M;
   _try(MatGetSize(A, &N, &M));
-
+ 
   // generalized eigenvalue problem A x - \lambda B x = 0
   EPS eps;
   _try(EPSCreate(PETSC_COMM_WORLD, &eps));
@@ -57,10 +57,12 @@ void eigenSolver::solve(int numEigenValues, std::string which)
   else if(which == "largest")
     _try(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE));
 
+
   // print info
   const EPSType type;
   _try(EPSGetType(eps, &type));
   Msg::Info("SLEPc solution method: %s", type);
+
   PetscInt nev;
   _try(EPSGetDimensions(eps, &nev, PETSC_NULL, PETSC_NULL));
   Msg::Info("SLEPc number of requested eigenvalues: %d", nev);