diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 74e694db024616c7775777c85a546fe276d5d504..c6ab9e37e7235568bb8397b657cd9a327f337de2 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -11,6 +11,7 @@
#include "Octree.h"
#include "gmshAssembler.h"
#include "gmshLaplace.h"
+#include "gmshCrossConf.h"
#include "gmshConvexCombination.h"
#include "gmshLinearSystemGmm.h"
#include "gmshLinearSystemCSR.h"
@@ -322,8 +323,12 @@ void GFaceCompound::parametrize() const
if(!oct){
coordinates.clear();
- parametrize(ITERU);
- parametrize(ITERV);
+
+ //parametrize(ITERU);
+ //parametrize(ITERV);
+
+ parametrize_conformal();
+
// checkOrientation();
computeNormals();
buildOct();
@@ -468,6 +473,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
gmshLinearSystemGmm<double> *_lsysb = new gmshLinearSystemGmm<double>;
//gmshLinearSystemCSRGmm<double> lsys;
_lsysb->setPrec(1.e-15);
+ _lsysb->setGmres(1);
if(Msg::GetVerbosity() == 99) _lsysb->setNoisy(2);
_lsys = _lsysb;
#else
@@ -784,7 +790,7 @@ void GFaceCompound::parametrize_conformal() const
{
MVertex *v1 = ordered[0];
- MVertex *v2 = ordered[1];
+ MVertex *v2 = ordered[ordered.size()/2];
myAssembler.fixVertex(v1, 0, 1, 0);
myAssembler.fixVertex(v1, 0, 2, 0);
myAssembler.fixVertex(v2, 0, 1, 1);
@@ -804,60 +810,20 @@ void GFaceCompound::parametrize_conformal() const
}
}
-#if 0
- Msg::Debug("Creating term %d dofs numbered %d fixed",
- myAssembler.sizeOfR(), myAssembler.sizeOfF());
- for(unsigned int i = 0; i < ordered.size(); i++){
- MVertex *v1 = ordered[i];
- MVertex *v2 = ordered[(i+1)%ordered.size()];
- myAssembler.assemble(v1,0,2,v1,0,1, 0.5);
- myAssembler.assemble(v1,0,2,v2,0,1, 0.5);
- myAssembler.assemble(v2,0,2,v1,0,1,-0.5);
- myAssembler.assemble(v2,0,2,v2,0,1,-0.5);
- myAssembler.assemble(v1,0,1,v1,0,2, 0.5);
- myAssembler.assemble(v1,0,1,v2,0,2, 0.5);
- myAssembler.assemble(v2,0,1,v1,0,2,-0.5);
- myAssembler.assemble(v2,0,1,v2,0,2,-0.5);
- }
- // internal BC's
- for(std::list<std::list<GEdge*> >::const_iterator itc = _interior_loops.begin();
- itc != _interior_loops.end(); ++itc){
- success = orderVertices(*itc, ordered, coords);
- if(!success){
- Msg::Error("Could not order vertices on boundary");
- return;
- }
- for(unsigned int i = 0; i < ordered.size(); i++){
- MVertex *v1 = ordered[i];
- MVertex *v2 = ordered[(i+1)%ordered.size()];
- myAssembler.assemble(v1,0,2,v1,0,1, 0.5);
- myAssembler.assemble(v1,0,2,v2,0,1, 0.5);
- myAssembler.assemble(v2,0,2,v1,0,1,-0.5);
- myAssembler.assemble(v2,0,2,v2,0,1,-0.5);
- myAssembler.assemble(v1,0,1,v1,0,2, 0.5);
- myAssembler.assemble(v1,0,1,v2,0,2, 0.5);
- myAssembler.assemble(v2,0,1,v1,0,2,-0.5);
- myAssembler.assemble(v2,0,1,v2,0,2,-0.5);
- }
- }
-#endif
-
- //gmshConvexCombinationTerm laplace(model(), &ONE, 1);
gmshFunction<double> ONE(1.0);
- gmshFunction<double> P5( .5);
- gmshFunction<double> M5(-0.5);
+ gmshFunction<double> MONE(-1.0 );
gmshLaplaceTerm laplace1(model(), &ONE, 1);
gmshLaplaceTerm laplace2(model(), &ONE, 2);
- gmshLaplaceTerm laplace12(model(), &M5, 1 , 2);
- gmshLaplaceTerm laplace21(model(), &P5, 2 , 1);
+ gmshCrossConfTerm cross12(model(), &MONE, 1 , 2);
+ gmshCrossConfTerm cross21(model(), &ONE, 2 , 1);
it = _compound.begin();
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
laplace1.addToMatrix(myAssembler, t);
laplace2.addToMatrix(myAssembler, t);
- laplace12.addToMatrix(myAssembler, t);
- laplace21.addToMatrix(myAssembler, t);
+ cross12.addToMatrix(myAssembler, t);
+ cross21.addToMatrix(myAssembler, t);
}
}
@@ -1145,7 +1111,7 @@ void GFaceCompound::getTriangle(double u, double v,
void GFaceCompound::buildOct() const
{
- // printStuff();
+ printStuff();
SBoundingBox3d bb;
int count = 0;
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index cc58bac166867bd809b02d6edcb7126ea345dee6..98786f31df43c2854942dcdc75173aef739f9336 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -7,6 +7,7 @@ set(SRC
Numeric.cpp
GmshMatrix.cpp
gmshLaplace.cpp
+ gmshCrossConf.cpp
gmshConvexCombination.cpp
gmshHelmholtz.cpp
gmshElasticity.cpp
diff --git a/Numeric/gmshCrossConf.cpp b/Numeric/gmshCrossConf.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6fd73a4340027ef82e87423fdf7d10d29b9aaf97
--- /dev/null
+++ b/Numeric/gmshCrossConf.cpp
@@ -0,0 +1,62 @@
+// 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>.
+
+#include "gmshCrossConf.h"
+#include "Numeric.h"
+
+void gmshCrossConfTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+{
+ int nbNodes = e->getNumVertices();
+ int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
+ int npts;
+ IntPt *GP;
+ double jac[3][3];
+ double invjac[3][3];
+ SVector3 Grads [256];
+ double grads[256][3];
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+
+ m.set_all(0.);
+
+ for (int i = 0; i < npts; i++){
+ const double u = GP[i].pt[0];
+ const double v = GP[i].pt[1];
+ const double w = GP[i].pt[2];
+ const double weight = GP[i].weight;
+ const double detJ = e->getJacobian(u, v, w, jac);
+ SPoint3 p; e->pnt(u, v, w, p);
+ const double _diff = (*_diffusivity)(p.x(), p.y(), p.z());
+ inv3x3(jac, invjac) ;
+ e->getGradShapeFunctions(u, v, w, grads);
+ for (int j = 0; j < nbNodes; j++){
+ Grads[j] = SVector3(invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
+ invjac[0][2] * grads[j][2],
+ invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
+ invjac[1][2] * grads[j][2],
+ invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
+ invjac[2][2] * grads[j][2]);
+ }
+
+ SVector3 N (jac[2][0],jac[2][1],jac[2][2]);
+ //SVector3 N (jac[0][2],jac[1][2],jac[2][2]);
+
+ double norm = N[0]*N[0]+ N[1]*N[1]+ N[2]*N[2];
+ if (norm != 1.0){
+ printf("Normal N=%g %g %g NORM = %g\n", N[0], N[1], N[2], norm);
+ //exit(1);
+ }
+
+ for (int j = 0; j < nbNodes; j++){
+ for (int k = 0; k <= j; k++){
+ m(j, k) += dot(crossprod(Grads[j],Grads[k]),N) * weight * detJ * _diff;
+ }
+ }
+ }
+
+ for (int j = 0; j < nbNodes; j++)
+ for (int k = 0; k < j; k++)
+ m(k, j) = -1.0*m(j, k);
+}
+
diff --git a/Numeric/gmshCrossConf.h b/Numeric/gmshCrossConf.h
new file mode 100644
index 0000000000000000000000000000000000000000..28ba7cc173daaded4e8941e99b9f4f9d386a4479
--- /dev/null
+++ b/Numeric/gmshCrossConf.h
@@ -0,0 +1,54 @@
+// 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 _GMSH_CROSSCONF_H_
+#define _GMSH_CROSSCONF_H_
+
+#include "gmshTermOfFormulation.h"
+#include "gmshFunction.h"
+#include "Gmsh.h"
+#include "GModel.h"
+#include "MElement.h"
+#include "GmshMatrix.h"
+
+class gmshCrossConfTerm : public gmshNodalFemTerm<double> {
+ protected:
+ const gmshFunction<double> *_diffusivity;
+ const int _iField;
+ int _iFieldB ;
+ public:
+ virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
+ virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
+ void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
+ {
+ *vR = e->getVertex(iRow);
+ *iCompR = 0;
+ *iFieldR = _iField;
+ }
+ void getLocalDofC(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
+ {
+ *vR = e->getVertex(iRow);
+ *iCompR = 0;
+ *iFieldR = _iFieldB;
+ }
+ public:
+ gmshCrossConfTerm(GModel *gm, gmshFunction<double> *diffusivity, int iField = 0, int iFieldB=-1) :
+ gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField)
+ {
+ _iFieldB = (iFieldB==-1) ? _iField : iFieldB;
+ }
+ virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+};
+
+class gmshCrossConfTerm2DParametric : gmshCrossConfTerm {
+ GFace *_gf;
+ public:
+ gmshCrossConfTerm2DParametric(GFace *gf, gmshFunction<double> *diffusivity, int iField = 0) :
+ gmshCrossConfTerm(gf->model(), diffusivity, iField), _gf(gf) {}
+ virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+};
+
+
+#endif
diff --git a/Numeric/gmshLaplace.cpp b/Numeric/gmshLaplace.cpp
index 63550e4744a76699f252ce8c4e2fb80083e0edca..9495e77183ce8cac03e8a7ed8c7717f6f8fc36b7 100644
--- a/Numeric/gmshLaplace.cpp
+++ b/Numeric/gmshLaplace.cpp
@@ -54,6 +54,4 @@ void gmshLaplaceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
m(k, j) = m(j, k);
}
-void gmshLaplaceTerm2DParametric::elementMatrix(MElement *e, gmshMatrix<double> &m) const
-{
-}
+
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
index e0428c8705ef6758ab0730c30bbdad1137c716b1..fb88b0b61af902f75bf2ad243195251b0db9784e 100644
--- a/Numeric/gmshLaplace.h
+++ b/Numeric/gmshLaplace.h
@@ -42,13 +42,6 @@ class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
};
-class gmshLaplaceTerm2DParametric : gmshLaplaceTerm {
- GFace *_gf;
- public:
- gmshLaplaceTerm2DParametric(GFace *gf, gmshFunction<double> *diffusivity, int iField = 0) :
- gmshLaplaceTerm(gf->model(), diffusivity, iField), _gf(gf) {}
- virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
-};
#endif
diff --git a/benchmarks/2d/square.geo b/benchmarks/2d/square.geo
index ec29fb7bdc8f39fcfa8fd9ff2fd4d09415cde094..7d9878cb983e9bd708dba7ad22a3bd9d34cabf3d 100644
--- a/benchmarks/2d/square.geo
+++ b/benchmarks/2d/square.geo
@@ -11,3 +11,4 @@ Line(4) = {1, 2};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(10) = {5};
+Compound Surface(11)={10};
diff --git a/benchmarks/boolean/TORUS.geo b/benchmarks/boolean/TORUS.geo
index 71b3e601ea4665bd80e956c36eeaa642fe6263b4..16be03c4a4fb6b4fa029e37195dc3a6e74a9cda9 100644
--- a/benchmarks/boolean/TORUS.geo
+++ b/benchmarks/boolean/TORUS.geo
@@ -1,3 +1,3 @@
OCCShape("Torus",{0,0,0,0,0,1,10,9},"None");
-BoundingBox{0,0,0,10,10,10};
+//BoundingBox{0,0,0,10,10,10};