From 9c6ea9bdce0a8d29cf6340d466ae2587a4b80ea1 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 2 Apr 2010 09:49:34 +0000
Subject: [PATCH]
---
Common/LuaBindings.cpp | 2 ++
Geo/GFaceCompound.cpp | 3 ++
Geo/MElement.cpp | 20 ++++++++++++
Geo/MElement.h | 1 +
Mesh/highOrderSmoother.cpp | 33 ++++++++------------
Solver/CMakeLists.txt | 1 +
Solver/dofManager.h | 2 +-
Solver/helmholtzTerm.h | 2 ++
Solver/linearSystem.cpp | 41 +++++++++++++++++++++++++
Solver/linearSystem.h | 3 ++
Solver/multiscaleLaplace.cpp | 59 ++++++++++++++++++++++++++++++++++--
11 files changed, 142 insertions(+), 25 deletions(-)
create mode 100644 Solver/linearSystem.cpp
diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 4451e9bbc1..48a0b7a5dd 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -25,6 +25,7 @@
#include "Bindings.h"
#include "drawContext.h"
#include "GmshMessage.h"
+#include "linearSystem.h"
#if defined(HAVE_SOLVER)
#include "linearSystemCSR.h"
@@ -379,6 +380,7 @@ binding::binding()
fullMatrix<double>::registerBindings(this);
gmshOptions::registerBindings(this);
Msg::registerBindings(this);
+ linearSystem<double>::registerBindings(this);
#if defined(HAVE_SOLVER)
function::registerBindings(this);
linearSystemCSRGmm<double>::registerBindings(this);
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 3e61adc6e4..14d50896f1 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -349,6 +349,8 @@ bool GFaceCompound::trivial() const
if(_compound.size() == 1 &&
(*(_compound.begin()))->getNativeType() == GEntity::OpenCascadeModel &&
(*(_compound.begin()))->geomType() != GEntity::DiscreteSurface){
+ if ((*(_compound.begin()))->periodic(0) ||
+ (*(_compound.begin()))->periodic(1) )return false;
return true;
}
return false;
@@ -798,6 +800,7 @@ 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
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index e58e1fe5be..8d055784ee 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -258,6 +258,26 @@ double MElement::getJacobian(double u, double v, double w, double jac[3][3])
return _computeDeterminantAndRegularize(this, jac);
}
+
+double MElement::getJacobian(double gsf[][3], double jac[3][3])
+{
+ jac[0][0] = jac[0][1] = jac[0][2] = 0.;
+ jac[1][0] = jac[1][1] = jac[1][2] = 0.;
+ jac[2][0] = jac[2][1] = jac[2][2] = 0.;
+
+ for (int i = 0; i < getNumVertices(); i++) {
+ const MVertex* v = getVertex(i);
+ double* gg = gsf[i];
+ for (int j = 0; j < 3; j++) {
+ jac[j][0] += v->x() * gg[j];
+ jac[j][1] += v->y() * gg[j];
+ jac[j][2] += v->z() * gg[j];
+ }
+ }
+
+ return _computeDeterminantAndRegularize(this, jac);
+}
+
double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
{
jac[0][0] = jac[0][1] = jac[0][2] = 0.;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index e903ec9c6d..d13f59b206 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -217,6 +217,7 @@ class MElement
int order=-1);
// return the Jacobian of the element evaluated at point (u,v,w) in
// parametric coordinates
+ double getJacobian(double gsf[][3], double jac[3][3]);
double getJacobian(double u, double v, double w, double jac[3][3]);
double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index cf41201773..7bdc487374 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -377,13 +377,13 @@ void highOrderSmoother::optimize(GFace * gf,
smooth(gf, true);
- /*
+
for (int i=0;i<100;i++){
int nbSwap =
swapHighOrderTriangles(gf,edgeVertices,faceVertices,this);
printf("%d swaps\n",nbSwap);
}
- */
+
// smooth(gf,true);
// smooth(gf,true);
// smooth(gf,true);
@@ -1382,15 +1382,11 @@ static int swapHighOrderTriangles(GFace *gf,
}
- for (edgeContainer::iterator ecit = edgeVertices.begin(); ecit != edgeVertices.end(); ecit++) {
- if (((*ecit).first.first == common_vertices[0] &&
- (*ecit).first.second == common_vertices[1]) ||
- ((*ecit).first.first == common_vertices[1] &&
- (*ecit).first.second == common_vertices[0])) {
- //printf ("Hehe, gotcha !\n");
- edgeVertices.erase(ecit);
- }
- }
+ std::vector<std::pair<MVertex*,MVertex*> > toDelete;
+ MEdge _temp(common_vertices[0],common_vertices[1]);
+ std::pair<MVertex*, MVertex*> _temp2(_temp.getMinVertex(), _temp.getMaxVertex());
+ edgeVertices.erase(_temp2);
+
/*
for (faceContainer::iterator fcit = faceVertices.begin(); fcit != faceVertices.end(); fcit++) {
bool remove_this = true;
@@ -1460,16 +1456,11 @@ static int swapHighOrderTriangles(GFace *gf,
break;
}
-
- for (edgeContainer::iterator ecit = edgeVertices.begin(); ecit != edgeVertices.end(); ecit++) {
- if (((*ecit).first.first == common_vertices[0] &&
- (*ecit).first.second == common_vertices[1]) ||
- ((*ecit).first.first == common_vertices[1] &&
- (*ecit).first.second == common_vertices[0])) {
- //printf ("Hehe, gotcha !\n");
- edgeVertices.erase(ecit);
- }
- }
+ // delete the edge inside the edgeContainer
+ std::vector<std::pair<MVertex*,MVertex*> > toDelete;
+ MEdge _temp(common_vertices[0],common_vertices[1]);
+ std::pair<MVertex*, MVertex*> _temp2(_temp.getMinVertex(), _temp.getMaxVertex());
+ edgeVertices.erase(_temp2);
/*
for (faceContainer::iterator fcit = faceVertices.begin(); fcit != faceVertices.end(); fcit++) {
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index 7194342b31..fbeea2c22e 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -4,6 +4,7 @@
# bugs and problems to <gmsh@geuz.org>.
set(SRC
+ linearSystem.cpp
linearSystemCSR.cpp
groupOfElements.cpp
elasticityTerm.cpp
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 49fa8f527d..9df7b57f88 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -108,7 +108,7 @@ class dofManager{
_linearSystems.insert(std::make_pair("B", l2));
}
inline void fixDof(Dof key, const dataVec &value)
- {
+ {
fixed[key] = value;
}
inline void fixDof(long int ent, int type, const dataVec &value)
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index 7bda156be0..cd2516c68e 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -52,6 +52,8 @@ class helmholtzTerm : public femTerm<scalar> {
// compute integration rule
const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() :
2 * (e->getPolynomialOrder() - 1);
+ // const int integrationOrder = 2 * e->getPolynomialOrder() + 1;
+
int npts; IntPt *GP;
e->getIntegrationPoints(integrationOrder, &npts, &GP);
// get the number of nodes
diff --git a/Solver/linearSystem.cpp b/Solver/linearSystem.cpp
new file mode 100644
index 0000000000..fa6c3f239b
--- /dev/null
+++ b/Solver/linearSystem.cpp
@@ -0,0 +1,41 @@
+#include "GmshConfig.h"
+#include "linearSystemFull.h"
+#include "linearSystemCSR.h"
+#include "linearSystemGmm.h"
+#include "linearSystemPETSc.h"
+
+#include "Bindings.h"
+
+template<>
+void linearSystem<double>::registerBindings(binding *b){
+ classBinding *cb = b->addClass<linearSystem<double> >("linearSystem");
+ cb->setDescription("An abstraction of a linear system Ax = b.");
+ methodBinding *cm;
+ cm = cb->addMethod("systemSolve", &linearSystem<double>::systemSolve);
+ cm->setDescription("compute x = A^{-1}b");
+
+#ifdef HAVE_TAUCS
+ cb = b->addClass<linearSystemCSRTaucs<double> >("linearSystemCSRTaucs");
+ cb->setDescription("A linear system solver, based on TAUCS, and that is ok for SDP sparse matrices.");
+ cm = cb->setConstructor<linearSystemCSRTaucs<double> >();
+ cm->setDescription ("A new TAUCS<double> solver");
+ cm->setArgNames(NULL);
+ cb->setParentClass<linearSystem<double> >();
+#endif
+
+#ifdef HAVE_PETSC
+ cb = b->addClass<linearSystemPETSc<double> >("linearSystemPETSc");
+ cb->setDescription("A linear system solver, based on PETSc");
+ cm = cb->setConstructor<linearSystemPETSc<double> >();
+ cm->setDescription ("A new PETSc<double> solver");
+ cm->setArgNames(NULL);
+ cb->setParentClass<linearSystem<double> >();
+#endif
+
+ cb = b->addClass<linearSystemFull<double> >("linearSystemFull");
+ cb->setDescription("A linear system solver, based on LAPACK (full matrices)");
+ cm = cb->setConstructor<linearSystemFull<double> >();
+ cm->setDescription ("A new Lapack based <double> solver");
+ cm->setArgNames(NULL);
+ cb->setParentClass<linearSystem<double> >();
+}
diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index c3e860aa05..3b6f9982ed 100644
--- a/Solver/linearSystem.h
+++ b/Solver/linearSystem.h
@@ -8,6 +8,8 @@
// A class that encapsulates a linear system solver interface :
// building a sparse matrix, solving a linear system
+class binding;
+
template <class scalar>
class linearSystem {
public :
@@ -24,6 +26,7 @@ class linearSystem {
virtual void zeroMatrix() = 0;
virtual void zeroRightHandSide() = 0;
virtual int systemSolve() = 0;
+ static void registerBindings (binding*);
};
#endif
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 1ea041f29f..1da275c77e 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -627,7 +627,8 @@ static void connected_left_right (std::vector<MElement *> &left,
static void printLevel(const char* fn,
std::vector<MElement *> &elements,
std::map<MVertex*,SPoint2> *coordinates,
- double version)
+ double version,
+ double *dx = 0)
{
if(!CTX::instance()->mesh.saveAll) return;
@@ -648,7 +649,10 @@ static void printLevel(const char* fn,
for (; it != vs.end() ; ++it){
(*it)->setIndex(index++);
SPoint2 p = (coordinates) ? (*coordinates)[*it] : SPoint2(0,0);
- if (coordinates) fprintf(fp, "%d %g %g 0\n", (*it)->getIndex(), p.x(), p.y());
+ if (coordinates) {
+ if (dx)fprintf(fp, "%d %22.15E %22.15E 0\n", (*it)->getIndex(), dx[2]*(p.x()-dx[0]), dx[2]*(p.y()-dx[1]));
+ else fprintf(fp, "%d %22.15E %22.15E 0\n", (*it)->getIndex(), p.x(), p.y());
+ }
else fprintf(fp, "%d %g %g %g\n", (*it)->getIndex(),
(*it)->x(), (*it)->y(), (*it)->z());
}
@@ -663,6 +667,10 @@ static void printLevel(const char* fn,
fclose(fp);
}
//--------------------------------------------------------------
+
+
+
+
static double localSize(MElement *e, std::map<MVertex*,SPoint2> &solution){
SBoundingBox3d local;
@@ -691,6 +699,35 @@ static double localSize(MElement *e, std::map<MVertex*,SPoint2> &solution){
}
+
+
+//--------------------------------------------------------------
+static void printLevel_onlysmall(const char* fn,
+ std::vector<MElement *> &elements,
+ std::map<MVertex*,SPoint2> *coordinates,
+ double version,
+ double tolerance){
+ std::vector<MElement *> small;
+ double dx[3] = {0,0,0};
+ int COUNT = 0;
+ for (int i=0;i<elements.size();i++){
+ double local_size = localSize(elements[i],*coordinates);
+ if (local_size < tolerance){
+ small.push_back(elements[i]);
+ for (int j=0;j<3;j++){
+ SPoint2 p = (*coordinates)[elements[i]->getVertex(j)];
+ dx[0] += p.x();
+ dx[1] += p.y();
+ COUNT++;
+ }
+ }
+ }
+ dx[0] /= COUNT;
+ dx[1] /= COUNT;
+ dx[2] = 1./tolerance;
+ printLevel(fn,small,coordinates,version,dx);
+}
+
//-------------------------------------------------------------
static void one2OneMap(std::vector<MElement *> &elements, std::map<MVertex*,SPoint2> &solution) {
@@ -790,6 +827,7 @@ 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);
@@ -947,7 +985,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
std::vector<SPoint2> localCoord;
double local_size = localSize(e,solution);
if (local_size < 1.e-5*global_size) //1.e-5
- tooSmall.push_back(e);
+ tooSmall.push_back(e);
else goodSize.push_back(e);
}
@@ -1068,8 +1106,10 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
//Save multiscale meshes
std::string name1(level._name+"real.msh");
std::string name2(level._name+"param.msh");
+ std::string name3(level._name+"param_small.msh");
printLevel (name1.c_str(),level.elements,0,2.2);
printLevel (name2.c_str(),level.elements,&level.coordinates,2.2);
+ printLevel_onlysmall (name3.c_str(),level.elements,&level.coordinates,2.2,1.e-15);
//For every small region compute a new parametrization
Msg::Info("Level (%d-%d): %d connected small regions",level.recur, level.region, regions.size());
@@ -1181,6 +1221,19 @@ void multiscaleLaplace::cut (std::vector<MElement *> &elements)
printLevel ("Rootcut-right.msh",right,0,2.2);
printLevel ("Rootcut-all.msh",elements, 0,2.2);
+ printLevel ("Rootcut-left-param.msh",left,&root->coordinates,2.2);
+ printLevel_onlysmall ("Rootcut-left-param10.msh",left,&root->coordinates,2.2,1.e-10);
+ printLevel_onlysmall ("Rootcut-left-param12.msh",left,&root->coordinates,2.2,1.e-12);
+ printLevel_onlysmall ("Rootcut-left-param15.msh",left,&root->coordinates,2.2,1.e-15);
+
+ printLevel ("Rootcut-right-param.msh",right,&root->coordinates,2.2);
+ printLevel_onlysmall ("Rootcut-right-param10.msh",right,&root->coordinates,2.2,1.e-10);
+ printLevel_onlysmall ("Rootcut-right-param12.msh",right,&root->coordinates,2.2,1.e-12);
+ printLevel_onlysmall ("Rootcut-right-param15.msh",right,&root->coordinates,2.2,1.e-15);
+
+ printLevel_onlysmall ("Rootcut-all-param12.msh",elements,&root->coordinates,2.2,1.e-12);
+ printLevel_onlysmall ("Rootcut-all-param15.msh",elements,&root->coordinates,2.2,1.e-15);
+
if ( elements.size() != left.size()+right.size()) {
Msg::Error("Cutting laplace wrong nb elements (%d) != left + right (%d)", elements.size(), left.size()+right.size());
exit(1);
--
GitLab