diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 20a69364cfdc3f33c733f4e6207014e5558f5c2a..1b9375be288ec841404251b409abd790c97f5b59 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -11,6 +11,7 @@
#include "Octree.h"
#include "gmshAssembler.h"
#include "gmshLaplace.h"
+#include "gmshConvexCombination.h"
#include "gmshLinearSystemGmm.h"
#include "gmshLinearSystemFull.h"
#include "FunctionSpace.h"
@@ -50,12 +51,14 @@ public:
class gmshDistanceBasedDiffusivity : public gmshFunction<double>
{
private:
+ mutable int comp;
MElement *_current;
mutable std::map<MVertex*, SPoint3> _coordinates;
public:
gmshDistanceBasedDiffusivity (std::map<MVertex*,SPoint3> &coordinates)
- : _current (0),_coordinates(coordinates){}
+ : _current (0),_coordinates(coordinates), comp(0){}
void setCurrent (MElement *current){ _current = current; }
+ void setCompound(int compound){ comp = compound; }
virtual double operator () (double x, double y, double z) const
{
double xyz[3] = {x, y, z}, uvw[3];
@@ -66,8 +69,11 @@ public:
value[i] = p[2];
}
double val = _current->interpolate(value, uvw[0], uvw[1], uvw[2]);
- //return exp(5*val);
+ //return 1./(exp(x)+1.e-4);//exp(5*val);
return 1.0;
+ //printf("compiound =%d \n", comp);
+ //if (comp == 8) return 1.0;
+ //else return 1.e-15;
}
};
@@ -98,13 +104,9 @@ void GFaceCompound::getBoundingEdges()
for (std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
- printf("set compound %d for face %d \n", tag(),(*it)->tag());
(*it)->setCompound(this);
}
- printf("***** In GFaceCompound: size U0=%d, v0=%d\n ", _U0.size(), _V0.size());
-
-
//in case the bounding edges are explicitely given
if (_U0.size()){
std::list<GEdge*> :: const_iterator it = _U0.begin();
@@ -146,7 +148,7 @@ void GFaceCompound::getBoundingEdges()
std::set<GEdge*>::iterator itf = _unique.begin();
for ( ; itf != _unique.end(); ++itf){
l_edges.push_back(*itf);
- printf("for edge %d, add face %d \n", (*itf)->tag(), this->tag());
+ //printf("for edge %d, add face %d \n", (*itf)->tag(), this->tag());
(*itf)->addFace(this);
}
@@ -491,12 +493,15 @@ void GFaceCompound::parametrize(iterationStep step) const
}
}
else{
+ //gmshConvexCombinationTerm laplace(model(), &ONE, 1);
gmshLaplaceTerm laplace(model(), &ONE, 1);
+ //gmshLaplaceTerm laplace(model(), &diffusivity, 1);
it = _compound.begin();
for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
- diffusivity.setCurrent(t);
+ //diffusivity.setCompound((*it)->tag());
+ //diffusivity.setCurrent(t);
laplace.addToMatrix(myAssembler, t);
}
}
@@ -917,7 +922,8 @@ void GFaceCompound::printStuff() const
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(),
- K1, K2, K3);
+ it0->second.z(),it1->second.z(),it2->second.z());
+ //K1, K2, K3);
fprintf(uvx,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
it0->second.x(), it0->second.y(), 0.0,
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index deb2a08c06b23539e4d560071e65c8ae560d2f37..525fc2edf89bbd77bdaf944db66a2db2fd7ee41d 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -916,7 +916,7 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
void GModel::createTopologyFromMesh()
{
- printf("***** In createTopologyFromMesh: \n");
+ Msg::Debug("***** In createTopologyFromMesh:");
std::vector<GEntity*> entities;
getEntities(entities);
@@ -944,10 +944,10 @@ void GModel::createTopologyFromMesh()
}
}
- printf("vertices size =%d \n", Dvertices.size());
- printf("edges size =%d \n", Dedges.size());
- printf("faces size =%d \n", Dfaces.size());
- printf("regions size =%d \n", Dregions.size());
+ Msg::Debug("vertices size =%d", Dvertices.size());
+ Msg::Debug("edges size =%d ", Dedges.size());
+ Msg::Debug("faces size =%d ", Dfaces.size());
+ Msg::Debug("regions size =%d", Dregions.size());
//For each discreteRegion, create Topology
//---------------------------------------
@@ -1094,7 +1094,6 @@ void GModel::createTopologyFromMesh()
vE = v1;
i=-1;
}
-
if (it == myEdges.end()) break;
}
printf("end Edges \n");
@@ -1132,7 +1131,6 @@ void GModel::createTopologyFromMesh()
if (std::find(all_vertices.begin(), all_vertices.end(), v1) == all_vertices.end()) all_vertices.push_back(v1);
}
e->mesh_vertices.insert(e->mesh_vertices.begin(), all_vertices.begin(), all_vertices.end());
- printf("all vertice size =%d\n", all_vertices.size());
for (std::vector<int>::iterator itFace = tagFaces.begin(); itFace != tagFaces.end(); itFace++) {
GFace *dFace = getFaceByTag(abs(*itFace));
diff --git a/Geo/GRegionCompound.cpp b/Geo/GRegionCompound.cpp
index 607e6c382825523939796a83076971517cf81294..0e6798612d138ceafa8f377ec4cb708ee1e14873 100644
--- a/Geo/GRegionCompound.cpp
+++ b/Geo/GRegionCompound.cpp
@@ -15,7 +15,7 @@
GRegionCompound::GRegionCompound(GModel *m, int tag, std::vector<GRegion*> &compound)
: GRegion(m, tag), _compound(compound)
{
- printf("********* In GRegion compound size =%d \n", _compound.size());
+ //printf("********* In GRegion compound size =%d \n", _compound.size());
getBoundingFaces();
}
@@ -29,7 +29,7 @@ void GRegionCompound::getBoundingFaces(){
std::multiset<GFace*> _touched;
std::vector<GRegion*>::iterator it = _compound.begin();
for ( ; it != _compound.end(); ++it){
- printf("face %d \n", (*it)->tag());
+ //printf("face %d \n", (*it)->tag());
std::list<GFace*> ed = (*it)->faces();
std::list<GFace*> :: iterator ite = ed.begin();
for ( ; ite != ed.end(); ++ite){
@@ -49,7 +49,7 @@ void GRegionCompound::getBoundingFaces(){
std::set<GFace*>::iterator itf = _unique.begin();
for ( ; itf != _unique.end(); ++itf){
l_faces.push_back(*itf);
- printf("for face %d, add region %d \n", (*itf)->tag(), this->tag());
+ //printf("for face %d, add region %d \n", (*itf)->tag(), this->tag());
(*itf)->addRegion(this);
}
diff --git a/Geo/Makefile b/Geo/Makefile
index a612fd0d282b2b8325fba0a4b390dc2ae1d5771b..c0c0a6cdbd0ecb050014bc5dbcfce73fb6659740 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -137,6 +137,8 @@ GFaceCompound${OBJEXT}: GFaceCompound.cpp ../Common/GmshConfig.h GFaceCompound.h
../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
../Geo/MElement.h ../Numeric/gmshFunction.h ../Common/Gmsh.h \
../Common/GmshMessage.h ../Numeric/GmshMatrix.h \
+ ../Numeric/gmshConvexCombination.h ../Numeric/gmshTermOfFormulation.h \
+ ../Numeric/gmshFunction.h ../Numeric/GmshMatrix.h \
../Numeric/gmshLinearSystemGmm.h ../Numeric/gmshLinearSystem.h \
../Numeric/gmshLinearSystemFull.h ../Numeric/gmshLinearSystem.h \
../Numeric/GmshMatrix.h
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 4eb5cc5c7d95ea187fc4cc3c84933ace94612ec0..22e518a7d21a5f2059eee6b9bf66745f81031c66 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -338,15 +338,6 @@ void discreteEdge::parametrize()
}
-
-
-// for (int i = 0; i < mesh_vertices.size(); i++){
-// double t1;
-// mesh_vertices[i]->getParameter(0,t1);
-// printf("** AFTER v1=%d t1=%g\n", mesh_vertices[i]->getNum(),t1 );
-// }
-
-
computeNormals();
}
@@ -387,8 +378,6 @@ void discreteEdge::computeNormals () const
//printf("normal =%g %g %g \n", itn->second.x(), itn->second.y(), itn->second.z());
}
- printf("********* normal size = %d \n", _normals.size());
-
}
void discreteEdge::getLocalParameter(const double &t, int &iLine,
diff --git a/Numeric/Makefile b/Numeric/Makefile
index 22a088a2dd38aa5ebee1f432d5da4dbb314ccfe9..d79c3d228e968d7e8492254b01aac8471aa6a5b0 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -14,6 +14,7 @@ CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE}
SRC = Numeric.cpp\
GmshMatrix.cpp\
gmshLaplace.cpp\
+ gmshConvexCombination.cpp\
gmshHelmholtz.cpp\
gmshElasticity.cpp\
gmshProjection.cpp\
@@ -79,6 +80,24 @@ gmshLaplace${OBJEXT}: gmshLaplace.cpp gmshLaplace.h gmshTermOfFormulation.h \
../Geo/MEdge.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
../Numeric/GmshMatrix.h ../Numeric/Gauss.h ../Common/Gmsh.h \
../Common/GmshMessage.h Numeric.h
+gmshConvexCombination${OBJEXT}: gmshConvexCombination.cpp \
+ gmshConvexCombination.h gmshTermOfFormulation.h GmshMatrix.h \
+ ../Common/GmshConfig.h ../Common/GmshMessage.h gmshFunction.h \
+ gmshAssembler.h gmshLinearSystem.h ../Geo/GModel.h ../Geo/GVertex.h \
+ ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+ ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/SOrientedBoundingBox.h \
+ ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/Pair.h \
+ ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
+ ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
+ ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
+ ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
+ ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+ ../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/MEdge.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
+ ../Numeric/GmshMatrix.h ../Numeric/Gauss.h ../Common/Gmsh.h \
+ ../Common/GmshMessage.h Numeric.h
gmshHelmholtz${OBJEXT}: gmshHelmholtz.cpp gmshHelmholtz.h \
gmshTermOfFormulation.h GmshMatrix.h ../Common/GmshConfig.h \
../Common/GmshMessage.h gmshFunction.h gmshAssembler.h \
diff --git a/Numeric/gmshConvexCombination.cpp b/Numeric/gmshConvexCombination.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..062755fc6fcc590d561957800b3c5ad66f6c3fe5
--- /dev/null
+++ b/Numeric/gmshConvexCombination.cpp
@@ -0,0 +1,28 @@
+// 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 "gmshConvexCombination.h"
+#include "Numeric.h"
+
+void gmshConvexCombinationTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+{
+
+ m.set_all(0.);
+ int nbNodes = e->getNumVertices();
+ double x = 0.3*(e->getVertex(0)->x()+e->getVertex(1)->x()+e->getVertex(2)->x());
+ double y = 0.3*(e->getVertex(0)->y()+e->getVertex(1)->y()+e->getVertex(2)->y());
+ double z = 0.3*(e->getVertex(0)->z()+e->getVertex(1)->z()+e->getVertex(2)->z());
+
+ const double _diff = 1.0;
+ //const double _diff = (*_diffusivity)(x, y, z);
+
+ for (int j = 0; j < nbNodes; j++){
+ for (int k = 0; k < nbNodes; k++){
+ m(j,k) = -1.*_diff;
+ }
+ m(j,j) = (nbNodes-1)*_diff;
+ }
+
+}
diff --git a/Numeric/gmshConvexCombination.h b/Numeric/gmshConvexCombination.h
new file mode 100644
index 0000000000000000000000000000000000000000..a2ae7c321bb1bde1924461bfce3186aaa838da64
--- /dev/null
+++ b/Numeric/gmshConvexCombination.h
@@ -0,0 +1,36 @@
+// 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_CONVEX_H_
+#define _GMSH_CONVEX_H_
+
+#include "gmshTermOfFormulation.h"
+#include "gmshFunction.h"
+#include "Gmsh.h"
+#include "GModel.h"
+#include "MElement.h"
+#include "GmshMatrix.h"
+
+class gmshConvexCombinationTerm : public gmshNodalFemTerm<double> {
+ protected:
+ const gmshFunction<double> *_diffusivity;
+ const int _iField ;
+ 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;
+ }
+ public:
+ gmshConvexCombinationTerm(GModel *gm, gmshFunction<double> *diffusivity, int iField = 0) :
+ gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField){}
+ virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+};
+
+
+#endif
diff --git a/Numeric/gmshLaplace.cpp b/Numeric/gmshLaplace.cpp
index 634ac96cb9dcb092f060065c6789442d2520c274..3eea76f3a8c6060e2efde2ffc8373311181a12de 100644
--- a/Numeric/gmshLaplace.cpp
+++ b/Numeric/gmshLaplace.cpp
@@ -38,14 +38,14 @@ void gmshLaplaceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
invjac[2][2] * grads[j][2];
}
double pi = 3.14;
- double H_x=1.0;
- double H_y=1.0;
- double H_z=1.0;
+ double K_x=1.0;
+ double K_y=1.0;
+ double K_z=1.0;
for (int j = 0; j < nbNodes; j++){
for (int k = 0; k <= j; k++){
- m(j, k) += (H_x*Grads[j][0] * Grads[k][0] +
- H_y*Grads[j][1] * Grads[k][1] +
- H_z*Grads[j][2] * Grads[k][2]) * weight * detJ * _diff;
+ m(j, k) += (K_x*Grads[j][0] * Grads[k][0] +
+ K_y*Grads[j][1] * Grads[k][1] +
+ K_z*Grads[j][2] * Grads[k][2]) * weight * detJ * _diff;
}
}
}
@@ -53,7 +53,27 @@ void gmshLaplaceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
for (int j = 0; j < nbNodes; j++)
for (int k = 0; k < j; k++)
m(k, j) = m(j, k);
-}
+
+
+ //check for positive scheme
+// //printf("ELEM\n");
+ for (int j = 0; j < nbNodes; j++){
+ double sum = m(j,0)+m(j,1)+m(j,2);
+ if (sum < 0.0){
+ printf("Sum LINE gmshLaplace Term NEG <0 %g %g %g\n", m(j,0),m(j,1),m(j,2));
+// for (int k = 0; k < nbNodes; k++) m(j,k) = -1.;
+// m(j,j) = (nbNodes-1);
+ }
+// else{
+// printf("Sum POS\n");
+// }
+// //for (int k = 0; k < nbNodes; k++) {
+// //printf("m(%d,%d)=%g ", j,k, m(j,k));
+// //}
+// //printf("\n");
+ }
+
+}
void gmshLaplaceTerm2DParametric::elementMatrix(MElement *e, gmshMatrix<double> &m) const
{