diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1ca56c7f262b2904b42cbbd386c8a012c1a106e5..6dfdf19cc88ae05313c788b5ce0f0697fbeea780 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -652,14 +652,14 @@ if(ENABLE_OCC)
endif(WIN32)
set(OCC_LIBS_REQUIRED
# subset of DataExchange
- TKSTEP TKSTEP209 TKSTEPAttr TKSTEPBase TKIGES TKXSBase
+ TKSTEP.0 TKSTEP209.0 TKSTEPAttr.0 TKSTEPBase.0 TKIGES.0 TKXSBase.0
# ModelingAlgorithms
- TKOffset TKFeat TKFillet TKBool TKShHealing TKMesh TKHLR TKBO TKPrim
- TKTopAlgo TKGeomAlgo
+ TKOffset.0 TKFeat.0 TKFillet.0 TKBool.0 TKShHealing.0 TKMesh.0 TKHLR.0 TKBO.0 TKPrim.0
+ TKTopAlgo.0 TKGeomAlgo.0
# ModelingData
- TKBRep TKGeomBase TKG3d TKG2d
+ TKBRep.0 TKGeomBase.0 TKG3d.0 TKG2d.0
# FoundationClasses
- TKAdvTools TKMath TKernel)
+ TKAdvTools.0 TKMath.0 TKernel.0)
list(LENGTH OCC_LIBS_REQUIRED NUM_OCC_LIBS_REQUIRED)
set(OCC_LIBS)
foreach(OCC ${OCC_LIBS_REQUIRED})
@@ -673,6 +673,7 @@ if(ENABLE_OCC)
endforeach(OCC)
list(LENGTH OCC_LIBS NUM_OCC_LIBS)
if(NUM_OCC_LIBS EQUAL NUM_OCC_LIBS_REQUIRED)
+
find_path(OCC_INC "BRep_Tool.hxx" PATHS ENV CASROOT PATH_SUFFIXES inc
include opencascade)
if(OCC_INC)
diff --git a/Common/GmshMessage.cpp b/Common/GmshMessage.cpp
index a7c6357060fc7cc4dff2fec3eb4b1339ac13f4fa..6e9fd668b7727ddba94086879157035a676ee1bc 100644
--- a/Common/GmshMessage.cpp
+++ b/Common/GmshMessage.cpp
@@ -66,15 +66,18 @@ static int vsnprintf(char *str, size_t size, const char *fmt, va_list ap)
void Msg::Init(int argc, char **argv)
{
#if defined(HAVE_MPI)
+ printf("initializing mpi\n");
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &_commRank);
MPI_Comm_size(MPI_COMM_WORLD, &_commSize);
MPI_Errhandler_set(MPI_COMM_WORLD, MPI_ERRORS_RETURN);
#endif
#if defined(HAVE_PETSC)
+ printf("initializing petsc\n");
PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
#endif
#if defined(HAVE_SLEPC)
+ printf("initializing slepc\n");
SlepcInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
#endif
time_t now;
diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 212d0a6764c02c7c700ab077fd144da831572ebf..a2be2e31f04f3dfb69cf113c2acf4c896141b3e9 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -424,8 +424,10 @@ binding::binding()
function::registerBindings(this);
linearSystem<double>::registerBindings(this);
linearSystemCSRGmm<double>::registerBindings(this);
+#if defined(HAVE_LUA)
elasticitySolverRegisterBindings(this);
#endif
+#endif
#if defined(HAVE_POST)
PView::registerBindings(this);
PViewData::registerBindings(this);
diff --git a/Common/gmshpy.i b/Common/gmshpy.i
index b924241e2af669343902a72f0cb7ff687bc720c3..ce22c6191880057819bc1ca51433653487de48ca 100644
--- a/Common/gmshpy.i
+++ b/Common/gmshpy.i
@@ -9,8 +9,10 @@
#include "function.h"
#include "dofManager.h"
#include "linearSystem.h"
+ #include "linearSystemFull.h"
#include "linearSystemPETSc.h"
#include "linearSystemPETSc.cpp" // needed for the specialization
+ #include "linearSystemCSR.h"
#include "GEntity.h"
#include "GVertex.h"
#include "GEdge.h"
@@ -35,6 +37,11 @@
#include "DivideAndConquer.h"
#include "Gmsh.h"
#include "functionPython.h"
+ #include "Context.h"
+ #include "SVector3.h"
+ #include "SPoint3.h"
+ #include "SPoint2.h"
+ #include "GPoint.h"
class errorHandler: public GmshMessage {
void operator()(std::string level, std::string message){
//const char *color = colorDefault;
@@ -58,6 +65,7 @@
namespace std {
%template(IntVector) vector<int>;
%template(DoubleVector) vector<double, std::allocator<double> >;
+ %template(VectorOfDoubleVector) vector<vector<double, std::allocator<double> > >;
%template(StringVector) vector<std::string, std::allocator<std::string> >;
%template(GEntityVector) vector<GEntity*, std::allocator<GEntity*> >;
%template(GVertexVector) vector<GVertex*, std::allocator<GVertex*> >;
@@ -71,12 +79,15 @@ namespace std {
%template(GFaceList) list<GFace*, std::allocator<GFace*> >;
}
+%include "Context.h"
%include "fullMatrix.h"
%include "dofManager.h"
%include "GModel.h"
%include "function.h"
%include "linearSystem.h"
+%include "linearSystemFull.h"
%include "linearSystemPETSc.h"
+%include "linearSystemCSR.h"
%include "GEntity.h"
%include "GVertex.h"
%include "GEdge.h"
@@ -100,12 +111,19 @@ namespace std {
%include "PViewFactory.h"
%include "DivideAndConquer.h"
%include "Gmsh.h"
+%include "SVector3.h"
+%include "SPoint3.h"
+%include "SPoint2.h"
+%include "GPoint.h"
%template(dofManagerDouble) dofManager<double>;
%template(linearSystemDouble) linearSystem<double>;
%template(linearSystemFullMatrixDouble) linearSystem<fullMatrix<double> >;
+%template(linearSystemFullDouble) linearSystemFull<double> ;
%template(linearSystemPETScDouble) linearSystemPETSc<double>;
+%template(linearSystemTAUCSDouble) linearSystemCSRTaucs<double>;
%template(fullMatrixDouble) fullMatrix<double>;
+%template(fullVectorDouble) fullVector<double>;
%template(linearSystemPETScBlockDouble) linearSystemPETSc<fullMatrix<double> >;
%include "functionPython.h"
diff --git a/Fltk/openglWindow.cpp b/Fltk/openglWindow.cpp
index 9e020b0fc82299c7bf39adee90de8e295e5235e5..edde01642a6f479cf23e45476e8ba8a13e2e0495 100644
--- a/Fltk/openglWindow.cpp
+++ b/Fltk/openglWindow.cpp
@@ -15,6 +15,7 @@
#include "MElement.h"
#include "Numeric.h"
#include "FlGui.h"
+//#include <GL/glu.h>
#include "drawContext.h"
#include "Context.h"
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 370050ce70511dad623c083444df1d8dca804929..6e7ffa83fd629995dda0323811748dc7a42de53f 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1205,9 +1205,13 @@ void GFace::moveToValidRange(SPoint2 &pt) const
}
}
+
void GFace::addLayersOfQuads(int nLayers, GVertex *gv, double hmin, double ratio){
- /*
+
+ SVector3 ez (0,0,1);
std::list<GEdgeLoop>::iterator it = edgeLoops.begin();
+ FILE *f = fopen ("coucou.pos","w");
+ fprintf(f,"View \"\"{\n");
for (; it != edgeLoops.end(); ++it){
bool found = false;
// look if this edge loop has the GVertex as an endpoint
@@ -1217,16 +1221,49 @@ void GFace::addLayersOfQuads(int nLayers, GVertex *gv, double hmin, double ratio
}
// we found an edge loop with the GVertex that was specified
if (found){
+ // first build a list of edges in the parametric space
+ std::vector<std::pair<MVertex*,SPoint2> > contour;
for (GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); ++it2){
GEdge *ge = it2->ge;
- for (int i=0;i<ge->lines.size();i++){
- SPoint2 p[2];
- reparamMeshEdgeOnFace ( ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),gf,p[0],p[1]);
+ SPoint2 p[2];
+ if (it2->_sign == 1){
+ for (int i=0;i<ge->lines.size();i++){
+ reparamMeshEdgeOnFace ( ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),this,p[0],p[1]);
+ contour.push_back(std::make_pair(ge->lines[i]->getVertex(0),p[0]));
+ }
+ }
+ else {
+ for (int i=ge->lines.size()-1;i>=0;i--){
+ reparamMeshEdgeOnFace ( ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),this,p[0],p[1]);
+ contour.push_back(std::make_pair(ge->lines[i]->getVertex(1),p[1]));
+ }
}
}
+ double hlayer = hmin;
+ for (int j=0;j<nLayers;j++){
+ for (int i=0;i<contour.size();i++){
+ SPoint2 p0 = contour[(i+0) % contour.size()].second;
+ SPoint2 p1 = contour[(i+1) % contour.size()].second;
+ SPoint2 p2 = contour[(i+2) % contour.size()].second;
+ SVector3 p0p1 (p1.x()-p0.x(),p1.y()-p0.y(),0.0);
+ SVector3 p1p2 (p2.x()-p1.x(),p2.y()-p1.y(),0.0);
+ SVector3 n01 = crossprod(ez,p0p1);
+ SVector3 n12 = crossprod(ez,p1p2);
+ SVector3 n = (n01+n12)*-0.5;
+ n.normalize();
+ double u = p1.x() + n.x() * hmin;
+ double v = p1.y() + n.y() * hmin;
+ GPoint gp = point(SPoint2(u,v));
+ _additional_vertices.push_back(new MFaceVertex(gp.x(),gp.y(),gp.z(),this,u,v));
+ fprintf(f,"SP(%g, %g, 0){1};\n",gp.x(),gp.y());
+ }
+ hlayer *= ratio;
+ hmin += hlayer;
+ }
+ fprintf(f,"};\n");
+ fclose(f);
}
- }
- */
+ }
}
@@ -1251,7 +1288,7 @@ void GFace::registerBindings(binding *b)
mb->setDescription("return the list of edges bounding this surface");
mb = cb->addMethod("addLayersOfQuads", &GFace::addLayersOfQuads);
mb->setDescription("insert layers of quads");
- mb->setArgNames("nLayers","startingVertex","hmin", "ratio", NULL);
+ mb->setArgNames("nLayers","startingVertex", "hmin", "ratio", NULL);
/* mb = cb->addMethod("addPolygon", &GFace::addPolygon);
mb->setDescription("insert a polygon mesh element");
mb->setArgNames("polygon", NULL);*/
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 22359ab83d7a818b3684e9ecc262c76188ecedee..d7eb3495e84e1e80378d97644192f80bf9fa59c0 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -347,6 +347,37 @@ GEdge *OCCFactory::addNURBS(GModel *gm, GVertex *start, GVertex *end,
return 0;
}
+/*
+GEdge *OCCFactory::addBezierSurface(GModel *gm,
+ std::vector<GEdge *> & wires, // four edges indeed
+ std::vector<std::vector<double> > points)
+{
+
+ TColgp_Array2OfPnt ctrlPoints(1, nbControlPoints + 2);
+ int index = 1;
+ ctrlPoints.SetValue(index++, gp_Pnt(start->x(), start->y(), start->z()));
+ for (int i = 0; i < nbControlPoints; i++) {
+ gp_Pnt aP(points[i][0],points[i][1],points[i][2]);
+ ctrlPoints.SetValue(index++, aP);
+ }
+
+
+ BRepBuilderAPI_MakeFace aGenerator (aBezierSurface);
+ BRepBuilderAPI_MakeWire wire_maker;
+ for (unsigned j = 0; j < wires.size(); j++) {
+ GEdge *ge = wires[j];
+ OCCEdge *occe = dynamic_cast<OCCEdge*>(ge);
+ if (occe){
+ wire_maker.Add(occe->getTopoDS_Edge());
+ }
+ }
+ TopoDS_Wire myWire = wire_maker.Wire();
+ aGenerator.Add (myWire);
+ aGenerator.Build();
+ TopoDS_Shape aResult = aGenerator.Shape();
+ return gm->_occ_internals->addFaceToModel(gm, TopoDS::Face(aResult));
+}
+*/
GEntity *OCCFactory::revolve(GModel *gm, GEntity* base,
std::vector<double> p1,
std::vector<double> p2, double angle)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 6a7ab364f0a37aecde9f0a27a046c03fa5eb2721..3150b52e676f3ff9ff4a09f33b48d51af96b5952 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -146,6 +146,19 @@ SPoint3 MElement::barycenter()
return p;
}
+double MElement::getVolume(){
+ int npts;
+ IntPt *pts;
+ getIntegrationPoints(getDim()*(getPolynomialOrder()-1), &npts, &pts);
+ double vol;
+ for (int i=0;i<npts;i++){
+ vol += getJacobianDeterminant(pts[i].pt[0],pts[i].pt[1],pts[i].pt[2])
+ * pts[i].weight;
+ }
+ return vol;
+}
+
+
int MElement::getVolumeSign()
{
double v = getVolume();
@@ -403,7 +416,7 @@ double MElement::interpolate(double val[], double u, double v, double w, int str
return sum;
}
-void MElement::interpolateGrad(double val[], double u, double v, double w, double f[3],
+void MElement::interpolateGrad(double val[], double u, double v, double w, double f[],
int stride, double invjac[3][3], int order)
{
double dfdu[3] = {0., 0., 0.};
@@ -427,7 +440,7 @@ void MElement::interpolateGrad(double val[], double u, double v, double w, doubl
}
}
-void MElement::interpolateCurl(double val[], double u, double v, double w, double f[3],
+void MElement::interpolateCurl(double val[], double u, double v, double w, double f[],
int stride, int order)
{
double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 5bab37d7b9416f4e20144678b12a6fd83b29a497..0dd0b769be2887c184f1979aff11c2a58fdd99b0 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -200,7 +200,7 @@ class MElement
virtual void revert(){}
// get volume of element
- virtual double getVolume(){ return 0.; }
+ virtual double getVolume();
// return sign of volume (+1 or -1) for 3D elements (or 0 if element
// has zero volume)
@@ -267,9 +267,9 @@ class MElement
// divergence) at point (u,v,w) in parametric coordinates
double interpolate(double val[], double u, double v, double w, int stride=1,
int order=-1);
- void interpolateGrad(double val[], double u, double v, double w, double f[3],
+ void interpolateGrad(double val[], double u, double v, double w, double f[],
int stride=1, double invjac[3][3]=0, int order=-1);
- void interpolateCurl(double val[], double u, double v, double w, double f[3],
+ void interpolateCurl(double val[], double u, double v, double w, double f[],
int stride=3, int order=-1);
double interpolateDiv(double val[], double u, double v, double w, int stride=3,
int order=-1);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 7cafc430e62c4b3c645292cafc7bce66fca16f88..0b00477c58180d87cdd65ba46842e280ed78f70f 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -917,6 +917,11 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
gf->meshStatistics.nbTriangle,
gf->meshStatistics.nbGoodQuality);
+ gf->mesh_vertices.insert(gf->mesh_vertices.end(),
+ gf->_additional_vertices.begin(),
+ gf->_additional_vertices.end());
+ gf->_additional_vertices.clear();
+
return true;
}
@@ -1500,7 +1505,7 @@ void deMeshGFace::operator() (GFace *gf)
}
// for debugging, change value from -1 to -100;
-int debugSurface = -1;
+int debugSurface = -100;
void meshGFace::operator() (GFace *gf)
{
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index f0503c57a43cd7e1dbeac32415fea1de20d270b4..404a1aad05376813d460dec3a9ae5469404e0d96 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1900,9 +1900,8 @@ void recombineIntoQuads(GFace *gf,
}
// give it a try : add one quad layer on the
-void addOneLayerOnContour(GFace *gf, GVertex *gv)
-{
/*
+void addOneLayerOnContour(GFace *gf, GVertex *gv){
, int nbLayers, double hplus, double factor){
// for each vertex
std::map<MVertex*,std::vector<MVertex*> >layers;
@@ -1991,8 +1990,8 @@ void addOneLayerOnContour(GFace *gf, GVertex *gv)
gf->quadrangles = newQuads;
}
}
- */
}
+*/
void quadsToTriangles(GFace *gf, double minqual = -10000.)
{
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index 200e9c5ab538d1fe16eb159a7bb276335e9f3080..a94db547f49806dccf589ddbf5a748495e0e3bee 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -56,7 +56,6 @@ void buildListOfEdgeAngle(e2t_cont adj, std::vector<edge_angle> &edges_detected,
std::vector<edge_angle> &edges_lonly);
void laplaceSmoothing(GFace *gf, int niter=1);
void edgeSwappingLawson(GFace *gf);
-void addOneLayerOnContour(GFace *gf, GVertex *gv);
enum swapCriterion {SWCR_DEL, SWCR_QUAL, SWCR_NORM, SWCR_CLOSE};
enum splitCriterion {SPCR_CLOSE, SPCR_QUAL, SPCR_ALLWAYS};
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index ea0caf004031d9261df498a6829c20a288b0c82f..d7ea21b5230911cd7251665f1dbfa46b449e2354 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -46,6 +46,15 @@ class fullVector
inline const scalar * getDataPtr() const { return _data; }
inline scalar operator () (int i) const { return _data[i]; }
inline scalar & operator () (int i) { return _data[i]; }
+ // set
+ inline void set(int r, scalar v){
+ #ifdef _DEBUG
+ if (r >= _r || r < 0)
+ Msg::Fatal("invalid index to access fullVector : %i (size = %i)",
+ r, _r);
+ #endif
+ (*this)(r) = v;
+ }
// operations
inline scalar norm() const
@@ -81,7 +90,7 @@ class fullVector
}
inline void setAll(const scalar &m)
{
- for(int i = 0; i < _r; i++) _data[i] = m;
+ for(int i = 0; i < _r; i++) set(i,m);
}
inline void setAll(const fullVector<scalar> &m)
{
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 633f97a7fbb48d4fda637b6abc569ef2eb25cbbe..31cbea5fe5956271037b87a6fa97f74a408afc24 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -723,11 +723,11 @@ PView* elasticitySolver::buildStressesView (const std::string postFileName)
#endif
-#if defined (HAVE_LUA)
#include "Bindings.h"
void elasticitySolverRegisterBindings(binding *b)
{
+#if defined (HAVE_LUA)
classBinding *cb;
cb = b->addClass<elasticitySolver> ("elasticitySolver");
cb->setDescription("A class that enables to solve elasticity problems");
@@ -795,6 +795,6 @@ void elasticitySolverRegisterBindings(binding *b)
cm = cb->setConstructor<elasticitySolver,GModel*,int>();
cm->setDescription ("A new elasticitySolver. The parameter is the unknowns tag");
cm->setArgNames("model","tag",NULL);
+#endif
}
-#endif
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index d8c498bd155d3998d4cf2aa1ab9f7b1231fad3ef..57b212ff7e0f2c090d64fca57bb812a1287762ea 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -116,5 +116,4 @@ class elasticitySolver
class binding;
void elasticitySolverRegisterBindings(binding *b);
-
#endif
diff --git a/benchmarks/2d/conge.geo b/benchmarks/2d/conge.geo
index 41bfe3d0777d0a3d1db10176b6ff1fee6818f18c..4ff78396df084af8a73caac7ab2d3cfc850991b8 100644
--- a/benchmarks/2d/conge.geo
+++ b/benchmarks/2d/conge.geo
@@ -78,4 +78,4 @@ Physical Line(25) = {12,13,14,15,16,17,18,19,20,10};
Physical Surface(26) = {22,24};
Physical Line(27) = {10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 14, 13, 12, 11};
Physical Line(28) = {17, 16, 20, 19, 18, 15};
-Recombine Surface {24, 22};
+//Recombine Surface {24, 22};