diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp index 216f7948fdeeab0805db3764648c31c186160a6a..4707da41035bc78fd344d077815f64deede0d0dc 100644 --- a/Geo/GFaceCompound.cpp +++ b/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)); diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h index 4c585496a917520ca445de045562b97a94057aa5..1314e92ca89806d659aa6545274e400b473d393d 100644 --- a/Geo/GFaceCompound.h +++ b/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; diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp index 45b80b95f822f42fc78da7d1fda5150412168284..e15e9ba3867de80b5baad22e838bb49746ea878f 100644 --- a/Mesh/meshGFace.cpp +++ b/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); diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp index f7f8da48797bc6961f28dbd4cbb31abf67c3f8ee..a55383c3d041d5c1a92ea6b212f0f025d3e0642f 100644 --- a/Mesh/meshGFaceLloyd.cpp +++ b/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() / diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h index d27de599cc1621d6c37232cc40aa6d929a059608..092bbdfb05cd655d5078b223d8260806d5e3298e 100644 --- a/Solver/convexCombinationTerm.h +++ b/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; */ +/* } */ +/* } */ + }; diff --git a/Solver/diagBCTerm.h b/Solver/diagBCTerm.h new file mode 100644 index 0000000000000000000000000000000000000000..ce59034a8b006070e572288f9f499db9e50ca2e4 --- /dev/null +++ b/Solver/diagBCTerm.h @@ -0,0 +1,61 @@ + +// 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 diff --git a/Solver/dofManager.h b/Solver/dofManager.h index 1edcad24391ebedcac1b4194cfeca9e62da9a2c3..5e0f0605728a6cd05eba0a7b6698aba175f3f6fa 100644 --- a/Solver/dofManager.h +++ b/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 = diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp index da58477d3b349275e389fbe0d07e13af061de92f..f162a08f5c1b785bd0dbfe77da400d6aff97eca1 100644 --- a/Solver/eigenSolver.cpp +++ b/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);