diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index 36e1680bf8a1aa921aaa95510e3f43597667e5af..2c03d02748e89de81931c54e42754cf37c4e20be 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -7,7 +7,6 @@
#include <time.h>
#include "GmshConfig.h"
#include "GmshDefines.h"
-#include "GmshPredicates.h"
#include "GModel.h"
#include "GmshMessage.h"
#include "OpenFile.h"
@@ -18,6 +17,7 @@
#include "Generator.h"
#include "Field.h"
#include "Context.h"
+#include "robustPredicates.h"
#include "meshPartition.h"
#include "GmshDaemon.h"
@@ -49,7 +49,7 @@ int GmshInitialize(int argc, char **argv)
#endif
// Initialize robust predicates
- gmsh::exactinit();
+ robustPredicates::exactinit();
if(dummy) delete dummy;
return 1;
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 242039a1079c726ca9e609fcaa9cb21e6829330a..47770f1c32b2746e6892d18152a6aa90d7e550dd 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -13,7 +13,7 @@
#include "MQuadrangle.h"
#include "MElementCut.h"
#include "VertexArray.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include "Numeric.h"
#include "EigSolve.h"
#include "GaussLegendre1D.h"
@@ -351,8 +351,8 @@ void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
ym /= (double)ndata;
zm /= (double)ndata;
- gmshMatrix<double> U(ndata, na), V(na, na);
- gmshVector<double> sigma(na);
+ fullMatrix<double> U(ndata, na), V(na, na);
+ fullVector<double> sigma(na);
for(int i = 0; i < ndata; i++) {
U(i, 0) = points[i].x() - xm;
U(i, 1) = points[i].y() - ym;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index d0387baef3d6a09da501101130b67ce5f11f00c2..b53bf878083ee08a05cfbae44f95aaed024a4daf 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -4,6 +4,7 @@
// bugs and problems to <gmsh@geuz.org>.
#include "GmshConfig.h"
+#include "GmshDefines.h"
#include "GFaceCompound.h"
#include "MLine.h"
#include "MTriangle.h"
@@ -19,9 +20,8 @@
#include "gmshLinearSystemGmm.h"
#include "gmshLinearSystemCSR.h"
#include "gmshLinearSystemFull.h"
-#include "FunctionSpace.h"
-#include "GmshDefines.h"
-#include "GmshPredicates.h"
+#include "functionSpace.h"
+#include "robustPredicates.h"
#include "meshGFaceOptimize.h"
#include "MElementCut.h"
#include "GEntity.h"
@@ -44,7 +44,7 @@ static void computeCGKernelPolygon(std::map<MVertex*,SPoint3> &coordinates,
// place at CG KERNEL polygon
int nbPts = cavV.size();
- gmshMatrix<double> u(100,2);
+ fullMatrix<double> u(100,2);
int i = 0;
for(std::vector<MVertex*>::iterator it = cavV.begin(); it != cavV.end(); it++){
SPoint3 vsp = coordinates[*it];
@@ -205,7 +205,7 @@ bool GFaceCompound::checkOrientation() const
double p1[2] = {v1[0],v1[1]};
double p2[2] = {v2[0],v2[1]};
double p3[2] = {v3[0],v3[1]};
- a_new = gmsh::orient2d(p1, p2, p3);
+ a_new = robustPredicates::orient2d(p1, p2, p3);
if(iter == 0) a_old=a_new;
if(a_new*a_old < 0.){
oriented = false;
@@ -247,7 +247,7 @@ bool GFaceCompound::checkCavity(std::vector<MElement*> &vTri) const
double p2[2] = {v2[0],v2[1]};
double p3[2] = {v3[0],v3[1]};
//printf("p1=(%g %g) p2=(%g %g) p3=(%g %g)\n",v1[0],v1[1], v2[0],v2[1], v3[0],v3[1] );
- a_new = gmsh::orient2d(p1, p2, p3);
+ a_new = robustPredicates::orient2d(p1, p2, p3);
if(i == 0) a_old=a_new;
if(a_new*a_old < 0.) badCavity = true;
a_old = a_new;
@@ -303,7 +303,7 @@ void GFaceCompound::one2OneMap() const
// double p1[2] = {v1[0],v1[1]};
// double p2[2] = {v2[0],v2[1]};
// double p3[2] = {v3[0],v3[1]};
-// double a_new = gmsh::orient2d(p1, p2, p3);
+// double a_new = robustPredicates::orient2d(p1, p2, p3);
// MVertex *e1=t->getVertex(0); MVertex *e2=t->getVertex(1); MVertex *e3=t->getVertex(2);
// printf("Point(%d)={%g, %g, 0}; \n ",e1->getNum(), v1[0],v1[1] );
// printf("Point(%d)={%g, %g, 0}; \n ",e2->getNum(), v2[0],v2[1] );
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 17d15258b1f9bf4df5e483f6a827ad82253ddaf3..d3d7d90984050c2ace89bc5b682024002db3bd0b 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2949,7 +2949,7 @@ struct PointSurface{
Surface *s;
};
-static void projectPS(gmshVector<double> &x, gmshVector<double> &res, void *data)
+static void projectPS(fullVector<double> &x, fullVector<double> &res, void *data)
{
PointSurface *ps = (PointSurface*)data;
Vertex c = InterpolateSurface(ps->s, x(0), x(1), 0, 0);
@@ -2967,7 +2967,7 @@ static void projectPS(gmshVector<double> &x, gmshVector<double> &res, void *data
bool ProjectPointOnSurface(Surface *s, Vertex &p, double uv[2])
{
- gmshVector<double> x(2);
+ fullVector<double> x(2);
x(0) = uv[0];
x(1) = uv[1];
PointSurface ps = {&p, s};
@@ -3151,7 +3151,7 @@ struct CurveSurface{
Surface *s;
};
-static void intersectCS(gmshVector<double> &uvt, gmshVector<double> &res, void *data)
+static void intersectCS(fullVector<double> &uvt, fullVector<double> &res, void *data)
{
CurveSurface *cs = (CurveSurface*)data;
Vertex vs = InterpolateSurface(cs->s, uvt(0), uvt(1), 0, 0);
@@ -3175,7 +3175,7 @@ bool IntersectCurvesWithSurface(List_T *curve_ids, int surface_id, List_T *shape
Curve *c = FindCurve((int)curve_id);
if(c){
CurveSurface cs = {c, s};
- gmshVector<double> uvt(3);
+ fullVector<double> uvt(3);
uvt(0) = 0.5;
uvt(1) = 0.5;
uvt(2) = 0.5;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 69ee7515c946b8191fd4a146ee5b38a3f4654dcd..5231fff5e4da3a4aafa25458818525c8299eae12 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -110,7 +110,7 @@ double MElement::rhoShapeMeasure()
void MElement::getShapeFunctions(double u, double v, double w, double s[], int o)
{
- const gmshFunctionSpace* fs = getFunctionSpace(o);
+ const functionSpace* fs = getFunctionSpace(o);
if(fs) fs->f(u, v, w, s);
else Msg::Error("Function space not implemented for this type of element");
}
@@ -118,7 +118,7 @@ void MElement::getShapeFunctions(double u, double v, double w, double s[], int o
void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3],
int o)
{
- const gmshFunctionSpace* fs = getFunctionSpace(o);
+ const functionSpace* fs = getFunctionSpace(o);
if(fs) fs->df(u, v, w, s);
else Msg::Error("Function space not implemented for this type of element");
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index d7daec74470df49a46c9c14aba64ac7df0cb02d5..54a76fc919b08886d8929dab2e82d7f9fb4a741f 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -10,11 +10,11 @@
#include <algorithm>
#include <string>
#include "GmshDefines.h"
+#include "GmshMessage.h"
#include "MVertex.h"
#include "MEdge.h"
#include "MFace.h"
-#include "GmshMessage.h"
-#include "FunctionSpace.h"
+#include "functionSpace.h"
#include "Gauss.h"
class GFace;
@@ -167,7 +167,7 @@ class MElement
virtual std::string getInfoString();
// get the function space for the element
- virtual const gmshFunctionSpace* getFunctionSpace(int order=-1) const { return 0; }
+ virtual const functionSpace* getFunctionSpace(int order=-1) const { return 0; }
// return the interpolating nodal shape functions evaluated at point
// (u,v,w) in parametric coordinates (if order == -1, use the
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index 7277f0fb4c1672e6800277d1fa3d1adb9ceeda63..7027e5d9afb947190848298cb8aaf7718c5b6b0d 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -103,7 +103,7 @@ class MPolyhedron : public MElement {
return vol;
}
virtual int getVolumeSign() { return (getVolume() >= 0) ? 1 : -1; }
- virtual const gmshFunctionSpace* getFunctionSpace(int order=-1) const
+ virtual const functionSpace* getFunctionSpace(int order=-1) const
{
return _orig->getFunctionSpace(order);
}
@@ -207,7 +207,7 @@ class MPolygon : public MElement {
}
virtual int getType() const { return TYPE_POLYG; }
virtual int getTypeForMSH() const { return MSH_POLYG_; }
- virtual const gmshFunctionSpace* getFunctionSpace(int order=-1) const
+ virtual const functionSpace* getFunctionSpace(int order=-1) const
{
return _orig->getFunctionSpace(order);
}
@@ -255,7 +255,7 @@ class MTriangleBorder : public MTriangle {
if(_domains[1]) return _domains[1]->getParent();
return NULL;
}
- virtual const gmshFunctionSpace* getFunctionSpace(int order=-1) const
+ virtual const functionSpace* getFunctionSpace(int order=-1) const
{
return getParent()->getFunctionSpace(order);
}
@@ -295,7 +295,7 @@ class MLineBorder : public MLine {
if(_domains[1]) return _domains[1]->getParent();
return NULL;
}
- virtual const gmshFunctionSpace* getFunctionSpace(int order=-1) const
+ virtual const functionSpace* getFunctionSpace(int order=-1) const
{
return getParent()->getFunctionSpace(order);
}
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 96a32cd12bb84aee11f2f800f428a6824ca1d56c..16c1a7a6dbc5096b5d76048e72d4c073b497c4da 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -8,16 +8,16 @@
#include "Context.h"
#include "qualityMeasures.h"
-const gmshFunctionSpace* MLine::getFunctionSpace(int o) const
+const functionSpace* MLine::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
- case 1: return &gmshFunctionSpaces::find(MSH_LIN_2);
- case 2: return &gmshFunctionSpaces::find(MSH_LIN_3);
- case 3: return &gmshFunctionSpaces::find(MSH_LIN_4);
- case 4: return &gmshFunctionSpaces::find(MSH_LIN_5);
- case 5: return &gmshFunctionSpaces::find(MSH_LIN_6);
+ case 1: return &functionSpaces::find(MSH_LIN_2);
+ case 2: return &functionSpaces::find(MSH_LIN_3);
+ case 3: return &functionSpaces::find(MSH_LIN_4);
+ case 4: return &functionSpaces::find(MSH_LIN_5);
+ case 5: return &functionSpaces::find(MSH_LIN_6);
default: Msg::Error("Order %d line function space not implemented", order);
}
return 0;
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 8f4689f5619befa55121ef5787d2812631c5a8aa..265c4662ebed76f85722c37ef36aca0e8a387bc8 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -63,7 +63,7 @@ class MLine : public MElement {
{
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
}
- virtual const gmshFunctionSpace* getFunctionSpace(int o=-1) const;
+ virtual const functionSpace* getFunctionSpace(int o=-1) const;
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index efc4d8d6eb2ea172a722d4d204e1c6e2f6f4d96e..f8b51640869a9eedce205cd31a8224a91b268f1d 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -60,7 +60,7 @@ void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3])
sys3x3(mat, b, uvw, &det);
}
-const gmshFunctionSpace* MTetrahedron::getFunctionSpace(int o) const
+const functionSpace* MTetrahedron::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
@@ -68,21 +68,21 @@ const gmshFunctionSpace* MTetrahedron::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 1: return &gmshFunctionSpaces::find(MSH_TET_4);
- case 2: return &gmshFunctionSpaces::find(MSH_TET_10);
- case 3: return &gmshFunctionSpaces::find(MSH_TET_20);
- case 4: return &gmshFunctionSpaces::find(MSH_TET_34);
- case 5: return &gmshFunctionSpaces::find(MSH_TET_52);
+ case 1: return &functionSpaces::find(MSH_TET_4);
+ case 2: return &functionSpaces::find(MSH_TET_10);
+ case 3: return &functionSpaces::find(MSH_TET_20);
+ case 4: return &functionSpaces::find(MSH_TET_34);
+ case 5: return &functionSpaces::find(MSH_TET_52);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return &gmshFunctionSpaces::find(MSH_TET_4);
- case 2: return &gmshFunctionSpaces::find(MSH_TET_10);
- case 3: return &gmshFunctionSpaces::find(MSH_TET_20);
- case 4: return &gmshFunctionSpaces::find(MSH_TET_35);
- case 5: return &gmshFunctionSpaces::find(MSH_TET_56);
+ case 1: return &functionSpaces::find(MSH_TET_4);
+ case 2: return &functionSpaces::find(MSH_TET_10);
+ case 3: return &functionSpaces::find(MSH_TET_20);
+ case 4: return &functionSpaces::find(MSH_TET_35);
+ case 5: return &functionSpaces::find(MSH_TET_56);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 1e13b294cb9e292c4dcc071fc45683598d3d5fa0..cf15ea0224882a3b4930dcffce604036fe58ecbc 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -129,7 +129,7 @@ class MTetrahedron : public MElement {
virtual double distoShapeMeasure();
virtual double etaShapeMeasure();
void xyz2uvw(double xyz[3], double uvw[3]);
- virtual const gmshFunctionSpace* getFunctionSpace(int o=-1) const;
+ virtual const functionSpace* getFunctionSpace(int o=-1) const;
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index a7331ce42905880eae83beade1f4a92b8dac24af..176d0725f2dafaabfb0d3d69d448e1e95bcbb2e1 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -30,7 +30,7 @@ double MTriangle::gammaShapeMeasure()
return qmTriangle(this, QMTRI_RHO);
}
-const gmshFunctionSpace* MTriangle::getFunctionSpace(int o) const
+const functionSpace* MTriangle::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
@@ -38,21 +38,21 @@ const gmshFunctionSpace* MTriangle::getFunctionSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 1: return &gmshFunctionSpaces::find(MSH_TRI_3);
- case 2: return &gmshFunctionSpaces::find(MSH_TRI_6);
- case 3: return &gmshFunctionSpaces::find(MSH_TRI_9);
- case 4: return &gmshFunctionSpaces::find(MSH_TRI_12);
- case 5: return &gmshFunctionSpaces::find(MSH_TRI_15I);
+ case 1: return &functionSpaces::find(MSH_TRI_3);
+ case 2: return &functionSpaces::find(MSH_TRI_6);
+ case 3: return &functionSpaces::find(MSH_TRI_9);
+ case 4: return &functionSpaces::find(MSH_TRI_12);
+ case 5: return &functionSpaces::find(MSH_TRI_15I);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return &gmshFunctionSpaces::find(MSH_TRI_3);
- case 2: return &gmshFunctionSpaces::find(MSH_TRI_6);
- case 3: return &gmshFunctionSpaces::find(MSH_TRI_10);
- case 4: return &gmshFunctionSpaces::find(MSH_TRI_15);
- case 5: return &gmshFunctionSpaces::find(MSH_TRI_21);
+ case 1: return &functionSpaces::find(MSH_TRI_3);
+ case 2: return &functionSpaces::find(MSH_TRI_6);
+ case 3: return &functionSpaces::find(MSH_TRI_10);
+ case 4: return &functionSpaces::find(MSH_TRI_15);
+ case 5: return &functionSpaces::find(MSH_TRI_21);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
}
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 13516f421cce53a15d59dc141f0d0c04541f3fa4..fef5d9d02e915856bbbd02e7d9db913e08e56d81 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -108,7 +108,7 @@ class MTriangle : public MElement {
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
- virtual const gmshFunctionSpace* getFunctionSpace(int o=-1) const;
+ virtual const functionSpace* getFunctionSpace(int o=-1) const;
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Geo/SOrientedBoundingBox.cpp b/Geo/SOrientedBoundingBox.cpp
index 838604f52ce8f7fc3fd88ee42fa663d19bab0c78..773240a51c70d48d64d8557dec25718550fafb1c 100644
--- a/Geo/SOrientedBoundingBox.cpp
+++ b/Geo/SOrientedBoundingBox.cpp
@@ -11,7 +11,7 @@
#include <time.h>
#include <stdlib.h>
#include "SOrientedBoundingBox.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include "DivideAndConquer.h"
#include "SBoundingBox3d.h"
@@ -198,7 +198,7 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
{
int num_vertices = vertices.size();
// First organize the data
- gmshMatrix<double> data(3,num_vertices);
+ fullMatrix<double> data(3,num_vertices);
int tmp = 0;
for (int i = 0; i < num_vertices; i++) {
bool okay = true;
@@ -220,7 +220,7 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
num_vertices = tmp;
// Compute the empirical means
- gmshVector<double> mean(3);
+ fullVector<double> mean(3);
for (int i = 0; i < 3; i++) {
mean(i) = 0;
for (int j = 0; j < num_vertices; j++) {
@@ -230,7 +230,7 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
}
// Get the deviation from the mean
- gmshMatrix<double> B(3,num_vertices);
+ fullMatrix<double> B(3,num_vertices);
for (int i = 0; i < 3; i++) {
for (int j = 0; j < num_vertices; j++) {
B(i,j) = data(i,j) - mean(i);
@@ -238,7 +238,7 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
}
// Compute the covariance matrix
- gmshMatrix<double> covariance(3,3);
+ fullMatrix<double> covariance(3,3);
B.mult(B.transpose(), covariance);
covariance.scale(1./(num_vertices-1));
/*
@@ -254,18 +254,18 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
}
}
- gmshMatrix<double> left_eigv(3,3);
- gmshMatrix<double> right_eigv(3,3);
- gmshVector<double> real_eig(3);
- gmshVector<double> img_eig(3);
+ fullMatrix<double> left_eigv(3,3);
+ fullMatrix<double> right_eigv(3,3);
+ fullVector<double> real_eig(3);
+ fullVector<double> img_eig(3);
covariance.eig(left_eigv, real_eig, img_eig, right_eigv,true);
// Now, project the data in the new basis.
- gmshMatrix<double> projected(3,num_vertices);
+ fullMatrix<double> projected(3,num_vertices);
left_eigv.transpose().mult(data,projected);
// Get the size of the box in the new direction
- gmshVector<double> mins(3);
- gmshVector<double> maxs(3);
+ fullVector<double> mins(3);
+ fullVector<double> maxs(3);
for (int i = 0; i < 3; i++) {
mins(i) = DBL_MAX;
maxs(i) = -DBL_MAX;
@@ -359,10 +359,10 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
// Now, examinate all the directions given by the edges of the convex hull
// to find the one that lets us build the least-area bounding rectangle for then
// points.
- gmshVector<double> axis(2);
+ fullVector<double> axis(2);
axis(0) = 1;
axis(1) = 0;
- gmshVector<double> axis2(2);
+ fullVector<double> axis2(2);
axis2(0) = 0;
axis2(1) = 1;
SOrientedBoundingRectangle least_rectangle;
@@ -373,7 +373,7 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
for (std::vector<Segment>::iterator seg = convex_hull.begin(); seg != convex_hull.end(); seg++) {
- gmshVector<double> segment(2);
+ fullVector<double> segment(2);
//segment(0) = record.points[(*seg).from].where.h - record.points[(*seg).to].where.h;
//segment(1) = record.points[(*seg).from].where.v - record.points[(*seg).to].where.v;
segment(0) = points[(*seg).from]->x() - points[(*seg).to]->x();
@@ -384,7 +384,7 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
double sine = axis(1)*segment(0) - segment(1)*axis(0);
//double sine = axis(0)*segment(1) - segment(0)*axis(1);
- gmshMatrix<double> rotation(2,2);
+ fullMatrix<double> rotation(2,2);
rotation(0,0) = cosine; rotation(0,1) = sine;
rotation(1,0) = -sine; rotation(1,1) = cosine;
@@ -395,12 +395,12 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
double min_y = DBL_MAX;
for (int i = 0; i < record.numPoints; i++) {
- gmshVector<double> pnt(2);
+ fullVector<double> pnt(2);
//pnt(0) = record.points[i].where.h;
//pnt(1) = record.points[i].where.v;
pnt(0) = points[i]->x();
pnt(1) = points[i]->y();
- gmshVector<double> rot_pnt(2);
+ fullVector<double> rot_pnt(2);
rotation.mult(pnt,rot_pnt);
if (rot_pnt(0) < min_x) min_x = rot_pnt(0);
if (rot_pnt(0) > max_x) max_x = rot_pnt(0);
@@ -409,19 +409,19 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
}
/**/
- gmshVector<double> center_rot(2);
- gmshVector<double> center_before_rot(2);
+ fullVector<double> center_rot(2);
+ fullVector<double> center_before_rot(2);
center_before_rot(0) = (max_x+min_x)/2.0;
center_before_rot(1) = (max_y+min_y)/2.0;
- gmshMatrix<double> rotation_inv(2,2);
+ fullMatrix<double> rotation_inv(2,2);
rotation_inv(0,0) = cosine; rotation_inv(0,1) = -sine;
rotation_inv(1,0) = sine; rotation_inv(1,1) = cosine;
rotation_inv.mult(center_before_rot,center_rot);
- gmshVector<double> axis_rot1(2);
- gmshVector<double> axis_rot2(2);
+ fullVector<double> axis_rot1(2);
+ fullVector<double> axis_rot2(2);
rotation_inv.mult(axis,axis_rot1);
rotation_inv.mult(axis2,axis_rot2);
diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index 8873bb7cade43b941dc47cd7209b0784dac0066c..979cf9e56262584c669907e14fb09387d45dce3f 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -12,8 +12,8 @@ void SMetric3::print (const char *s) const
SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
{
SMetric3 im1 = m1.invert();
- gmshMatrix<double> V(3,3);
- gmshVector<double> S(3);
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
im1 *= m2;
im1.eig(V,S,true);
SVector3 v0(V(0,0),V(1,0),V(2,0));
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index 0c3d0982c2b88eddbdd17aa4d7aa359c7f4154a3..af6c16471b08ddd21bf955e065065af756e0dcd3 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -7,7 +7,7 @@
#define _STENSOR3_H_
#include "SVector3.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include "Numeric.h"
// concrete class for symmetric positive definite 3x3 matrix
@@ -21,14 +21,14 @@ class SMetric3 {
static int _index[3][3] = {{0,1,3},{1,2,4},{3,4,5}};
return _index[i][j];
}
- void getMat (gmshMatrix<double> & mat) const{
+ void getMat (fullMatrix<double> & mat) const{
for (int i=0;i<3;i++){
for (int j=0;j<3;j++){
mat(i,j) = _val[getIndex(i,j)];
}
}
}
- void setMat (const gmshMatrix<double> & mat){
+ void setMat (const fullMatrix<double> & mat){
for (int i=0;i<3;i++)
for (int j=i;j<3;j++)
_val[getIndex(i,j)] = mat(i,j);
@@ -71,13 +71,13 @@ class SMetric3 {
// where the elements of diag are l_i = h_i^-2
// and the rows of e are the UNIT and ORTHOGONAL directions
- gmshMatrix<double> e(3,3);
+ fullMatrix<double> e(3,3);
e(0,0) = t1(0); e(0,1) = t1(1); e(0,2) = t1(2);
e(1,0) = t2(0); e(1,1) = t2(1); e(1,2) = t2(2);
e(2,0) = t3(0); e(2,1) = t3(1); e(2,2) = t3(2);
e.invertInPlace();
- gmshMatrix<double> tmp(3,3);
+ fullMatrix<double> tmp(3,3);
tmp(0,0) = l1 * e(0,0);
tmp(0,1) = l2 * e(0,1);
tmp(0,2) = l3 * e(0,2);
@@ -104,7 +104,7 @@ class SMetric3 {
return _val[getIndex(i,j)];
}
SMetric3 invert () const {
- gmshMatrix<double> m(3,3);
+ fullMatrix<double> m(3,3);
getMat(m);
m.invertInPlace();
SMetric3 ithis;
@@ -125,26 +125,26 @@ class SMetric3 {
return *this;
}
SMetric3& operator *= (const SMetric3 &other) {
- gmshMatrix<double> m1(3,3),m2(3,3),m3(3,3);
+ fullMatrix<double> m1(3,3),m2(3,3),m3(3,3);
getMat(m1);
other.getMat(m2);
m1.mult(m2,m3);
setMat(m3);
return *this;
}
- SMetric3 transform (gmshMatrix<double> &V){
- gmshMatrix<double> m(3,3);
+ SMetric3 transform (fullMatrix<double> &V){
+ fullMatrix<double> m(3,3);
getMat(m);
- gmshMatrix<double> result(3,3),temp(3,3);
+ fullMatrix<double> result(3,3),temp(3,3);
V.transpose().mult(m,temp);
temp.mult(V,result);
SMetric3 a; a.setMat(result);
return a;
}
// s: true if eigenvalues are sorted (from min to max of the REAL part)
- void eig (gmshMatrix<double> &V, gmshVector<double> &S, bool s=false) const {
- gmshMatrix<double> me(3,3),right(3,3);
- gmshVector<double> im(3);
+ void eig (fullMatrix<double> &V, fullVector<double> &S, bool s=false) const {
+ fullMatrix<double> me(3,3),right(3,3);
+ fullVector<double> im(3);
getMat(me);
me.eig(V,S,im,right,s);
}
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index a341dcd82732df7d24a1313a9d707e0aa0273447..232ffaeea7016fb1df85cb1a1cb8a22c7cc6e5ad 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -6,7 +6,7 @@
#include <math.h>
#include <stdio.h>
#include "GmshMessage.h"
-#include "GmshPredicates.h"
+#include "robustPredicates.h"
#include "Numeric.h"
#include "BDS.h"
#include "GFace.h"
@@ -210,10 +210,10 @@ int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
// double q1[2] = {x3,y3};
// double q2[2] = {x4,y4};
-// double rp1 = gmsh::orient2d(p1, p2, q1);
-// double rp2 = gmsh::orient2d(p1, p2, q2);
-// double rq1 = gmsh::orient2d(q1, q2, p1);
-// double rq2 = gmsh::orient2d(q1, q2, p2);
+// double rp1 = robustPredicates::orient2d(p1, p2, q1);
+// double rp2 = robustPredicates::orient2d(p1, p2, q2);
+// double rq1 = robustPredicates::orient2d(q1, q2, p1);
+// double rq2 = robustPredicates::orient2d(q1, q2, p2);
// if(rp1*rp2<=0 && rq1*rq2<=0){
// // printf("p1 %22.15E %22.15E %22.15E %22.15E \n",x1,y1,x2,y2);
@@ -731,8 +731,8 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2,
double op1[2] = {_q1->u,_q1->v};
double op2[2] = {_q2->u,_q2->v};
- double ori_t1 = gmsh::orient2d(op1, p1, op2);
- double ori_t2 = gmsh::orient2d(op1, op2, p2);
+ double ori_t1 = robustPredicates::orient2d(op1, p1, op2);
+ double ori_t2 = robustPredicates::orient2d(op1, op2, p2);
// printf("%d %d %d %d %g %g\n",_p1->iD,_p2->iD,_q1->iD,_q2->iD,ori_t1,ori_t2);
return (ori_t1 * ori_t2 > 0); // the quadrangle was strictly convex !
@@ -978,8 +978,8 @@ static bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BD
double pc[2] = {pts[2]->u, pts[2]->v};
double pd[2] = {pts[3]->u, pts[3]->v};
- double ori_init1 = gmsh::orient2d(pa, pb, pc);
- double ori_init2 = gmsh::orient2d(pc, pd, pa);
+ double ori_init1 = robustPredicates::orient2d(pa, pb, pc);
+ double ori_init2 = robustPredicates::orient2d(pc, pd, pa);
if(p == pts[0]){
pa[0] = u;
@@ -1002,8 +1002,8 @@ static bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BD
return false;
}
- double ori_final1 = gmsh::orient2d(pa, pb, pc);
- double ori_final2 = gmsh::orient2d(pc, pd, pa);
+ double ori_final1 = robustPredicates::orient2d(pa, pb, pc);
+ double ori_final2 = robustPredicates::orient2d(pc, pd, pa);
// allow to move a point when a triangle was flat
return (ori_init1 * ori_final1 > 0) && (ori_init2 * ori_final2 > 0);
}
@@ -1026,7 +1026,7 @@ static bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v
if(area_init == 0.0) return true;
- double ori_init = gmsh::orient2d(pa, pb, pc);
+ double ori_init = robustPredicates::orient2d(pa, pb, pc);
if(p == pts[0]){
pa[0] = u;
@@ -1048,7 +1048,7 @@ static bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v
double area_final = fabs(a[0] * b[1] - a[1] * b[0]);
if(area_final < 0.1 * area_init) return false;
- double ori_final = gmsh::orient2d(pa, pb, pc);
+ double ori_final = robustPredicates::orient2d(pa, pb, pc);
// allow to move a point when a triangle was flat
return ori_init * ori_final > 0;
}
diff --git a/Mesh/DivideAndConquer.cpp b/Mesh/DivideAndConquer.cpp
index 486502d92353fe5ce8f8a9cfb4685c99f9ed8095..73504e2784887a5795852fc29d9da986b93a6cc7 100644
--- a/Mesh/DivideAndConquer.cpp
+++ b/Mesh/DivideAndConquer.cpp
@@ -16,9 +16,9 @@
// value to avoid 3 aligned points or 4 cocyclical points!
#include "GmshMessage.h"
-#include "GmshPredicates.h"
-#include "Numeric.h"
#include "DivideAndConquer.h"
+#include "Numeric.h"
+#include "robustPredicates.h"
#include "MallocUtils.h"
#define Pred(x) ((x)->prev)
@@ -85,7 +85,7 @@ int DocRecord::IsLeftOf(PointNumero x, PointNumero y, PointNumero check)
double pc[2] = {(double)points[check].where.h, (double)points[check].where.v};
// we use robust predicates here
- double result = gmsh::orient2d(pa, pb, pc);
+ double result = robustPredicates::orient2d(pa, pb, pc);
return result > 0;
}
@@ -168,7 +168,8 @@ int DocRecord::Qtest(PointNumero h, PointNumero i, PointNumero j, PointNumero k)
double pd[2] = {(double)points[k].where.h, (double)points[k].where.v};
// we use robust predicates here
- double result = gmsh::incircle(pa, pb, pc, pd) * gmsh::orient2d(pa, pb, pc);
+ double result = robustPredicates::incircle(pa, pb, pc, pd) *
+ robustPredicates::orient2d(pa, pb, pc);
return (result < 0) ? 1 : 0;
}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 8fa38bb91ab7964789289e115ec258520a2c80c0..9f055f443cb5ac4e67006ef2bebdca1cfa159f9f 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -20,8 +20,8 @@
#include "OS.h"
#include "Numeric.h"
#include "Context.h"
-#include "GmshMatrix.h"
-#include "FunctionSpace.h"
+#include "fullMatrix.h"
+#include "functionSpace.h"
#define SQU(a) ((a)*(a))
@@ -46,7 +46,7 @@ static bool mappingIsInvertible(MTetrahedron *e)
if (det0 * detN <= 0.) return false;
}
- const gmshMatrix<double> &points = e->getFunctionSpace()->points;
+ const fullMatrix<double> &points = e->getFunctionSpace()->points;
for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
const double u = points(i,0);
@@ -68,7 +68,7 @@ static double mylength(GEdge *ge, int i, double *u)
return ge->length(u[i], u[i+1], 10);
}
-static void myresid(int N, GEdge *ge, double *u, gmshVector<double> &r)
+static void myresid(int N, GEdge *ge, double *u, fullVector<double> &r)
{
double L[100];
for (int i = 0; i < N - 1; i++) L[i] = mylength(ge, i, u);
@@ -95,10 +95,10 @@ static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
// create the tangent matrix
const int M = N - 2;
- gmshMatrix<double> J(M, M);
- gmshVector<double> DU(M);
- gmshVector<double> R(M);
- gmshVector<double> Rp(M);
+ fullMatrix<double> J(M, M);
+ fullVector<double> DU(M);
+ fullVector<double> R(M);
+ fullVector<double> Rp(M);
int iter = 1;
@@ -140,7 +140,7 @@ static double mylength(GFace *gf, int i, double *u, double *v)
return gf->length(SPoint2(u[i], v[i]), SPoint2(u[i + 1], v[i + 1]), 10);
}
-static void myresid(int N, GFace *gf, double *u, double *v, gmshVector<double> &r)
+static void myresid(int N, GFace *gf, double *u, double *v, fullVector<double> &r)
{
double L[100];
for (int i = 0; i < N - 1; i++) L[i] = mylength(gf, i, u, v);
@@ -174,10 +174,10 @@ static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
// create the tangent matrix
const int M = N - 2;
- gmshMatrix<double> J(M, M);
- gmshVector<double> DU(M);
- gmshVector<double> R(M);
- gmshVector<double> Rp(M);
+ fullMatrix<double> J(M, M);
+ fullVector<double> DU(M);
+ fullVector<double> R(M);
+ fullVector<double> Rp(M);
int iter = 1;
@@ -409,20 +409,20 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
gf->geomType() == GEntity::BoundaryLayerSurface)
linear = true;
- gmshMatrix<double> points;
+ fullMatrix<double> points;
int start = 0;
switch (nPts){
case 2 :
- points = gmshFunctionSpaces::find(MSH_TRI_10).points;
+ points = functionSpaces::find(MSH_TRI_10).points;
start = 9;
break;
case 3 :
- points = gmshFunctionSpaces::find(MSH_TRI_15).points;
+ points = functionSpaces::find(MSH_TRI_15).points;
start = 12;
break;
case 4 :
- points = gmshFunctionSpaces::find(MSH_TRI_21).points;
+ points = functionSpaces::find(MSH_TRI_21).points;
start = 15;
break;
default :
@@ -559,20 +559,20 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
faceContainer &faceVertices, edgeContainer &edgeVertices,
bool linear, int nPts = 1)
{
- gmshMatrix<double> points;
+ fullMatrix<double> points;
int start = 0;
switch (nPts){
case 2 :
- points = gmshFunctionSpaces::find(MSH_TRI_10).points;
+ points = functionSpaces::find(MSH_TRI_10).points;
start = 9;
break;
case 3 :
- points = gmshFunctionSpaces::find(MSH_TRI_15).points;
+ points = functionSpaces::find(MSH_TRI_15).points;
start = 12;
break;
case 4 :
- points = gmshFunctionSpaces::find(MSH_TRI_21).points;
+ points = functionSpaces::find(MSH_TRI_21).points;
start = 15;
break;
default :
@@ -652,16 +652,16 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
std::vector<MVertex*> &vr, bool linear, int nPts = 1)
{
- gmshMatrix<double> points;
+ fullMatrix<double> points;
int start = 0;
switch (nPts){
case 3 :
- points = gmshFunctionSpaces::find(MSH_TET_35).points;
+ points = functionSpaces::find(MSH_TET_35).points;
start = 34;
break;
case 4 :
- points = gmshFunctionSpaces::find(MSH_TET_56).points;
+ points = functionSpaces::find(MSH_TET_56).points;
start = 52;
break;
default :
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index 9748f815cb31785141932c54d9a5e5bb5401b288..f5987a4eb3189870691aae0f09cd4154d30e7cec 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -102,7 +102,7 @@ struct p2data{
GFace *gf;
MTriangle *t1,*t2;
MVertex *n12;
- gmshMatrix<double> *m1,*m2;
+ fullMatrix<double> *m1,*m2;
gmshHighOrderSmoother *s;
p2data(GFace *_gf, MTriangle *_t1, MTriangle *_t2, MVertex *_n12,
gmshHighOrderSmoother *_s)
@@ -111,9 +111,9 @@ struct p2data{
elasticityTerm el(0, 1.e3, .3333,1);
s->moveToStraightSidedLocation(t1);
s->moveToStraightSidedLocation(t2);
- m1 = new gmshMatrix<double>(3 * t1->getNumVertices(),
+ m1 = new fullMatrix<double>(3 * t1->getNumVertices(),
3 * t1->getNumVertices());
- m2 = new gmshMatrix<double>(3 * t2->getNumVertices(),
+ m2 = new fullMatrix<double>(3 * t2->getNumVertices(),
3 * t2->getNumVertices());
el.elementMatrix(t1,*m1);
el.elementMatrix(t2,*m2);
@@ -130,7 +130,7 @@ struct pNdata{
GFace *gf;
MTriangle *t1,*t2;
const std::vector<MVertex*> &n;
- gmshMatrix<double> *m1,*m2;
+ fullMatrix<double> *m1,*m2;
gmshHighOrderSmoother *s;
pNdata(GFace *_gf,MTriangle *_t1,MTriangle *_t2, const std::vector<MVertex*> &_n,
gmshHighOrderSmoother *_s)
@@ -139,9 +139,9 @@ struct pNdata{
elasticityTerm el(0, 1.e3, .3333,1);
s->moveToStraightSidedLocation(t1);
s->moveToStraightSidedLocation(t2);
- m1 = new gmshMatrix<double>(3 * t1->getNumVertices(),
+ m1 = new fullMatrix<double>(3 * t1->getNumVertices(),
3 * t1->getNumVertices());
- m2 = new gmshMatrix<double>(3 * t2->getNumVertices(),
+ m2 = new fullMatrix<double>(3 * t2->getNumVertices(),
3 * t2->getNumVertices());
el.elementMatrix(t1,*m1);
el.elementMatrix(t2,*m2);
@@ -155,11 +155,11 @@ struct pNdata{
};
static double _DeformationEnergy(MElement *e,
- gmshMatrix<double> *K,
+ fullMatrix<double> *K,
gmshHighOrderSmoother *s)
{
int N = e->getNumVertices();
- gmshVector<double> Kdx(N*3),dx(N*3);
+ fullVector<double> Kdx(N*3),dx(N*3);
dx.scale(0.0);
Kdx.scale(0.0);
@@ -194,7 +194,7 @@ static double _DeformationEnergy(pNdata *p)
_DeformationEnergy(p->t2, p->m2, p->s);
}
-static double _function_p2tB(gmshVector<double> &x, void *data)
+static double _function_p2tB(fullVector<double> &x, void *data)
{
p2data *p = (p2data*)data;
GPoint gp12 = p->gf->point(SPoint2(x(0), x(1)));
@@ -212,7 +212,7 @@ static double _function_p2tB(gmshVector<double> &x, void *data)
return E;
}
-static double _function_p2t(gmshVector<double> &x, void *data)
+static double _function_p2t(fullVector<double> &x, void *data)
{
p2data *p = (p2data*)data;
GPoint gp12 = p->gf->point(SPoint2(x(0), x(1)));
@@ -229,7 +229,7 @@ static double _function_p2t(gmshVector<double> &x, void *data)
return -q;
}
-static double _function_pNt(gmshVector<double> &x, void *data)
+static double _function_pNt(fullVector<double> &x, void *data)
{
pNdata *p = (pNdata*)data;
static double xx[256];
@@ -253,7 +253,7 @@ static double _function_pNt(gmshVector<double> &x, void *data)
return -q;
}
-static double _function_pNtB(gmshVector<double> &x, void *data)
+static double _function_pNtB(fullVector<double> &x, void *data)
{
pNdata *p = (pNdata*)data;
static double xx[256];
@@ -390,9 +390,9 @@ void gmshHighOrderSmoother::optimize(GFace * gf,
void gmshHighOrderSmoother::computeMetricVector(GFace *gf,
MElement *e,
- gmshMatrix<double> &J,
- gmshMatrix<double> &JT,
- gmshVector<double> &D)
+ fullMatrix<double> &J,
+ fullMatrix<double> &JT,
+ fullVector<double> &D)
{
int nbNodes = e->getNumVertices();
for (int j = 0; j < nbNodes; j++){
@@ -556,13 +556,13 @@ double gmshHighOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
const int n2 = 2 * nbNodes;
const int n3 = 3 * nbNodes;
- gmshMatrix<double> K33(n3, n3);
- gmshMatrix<double> K22(n2, n2);
- gmshMatrix<double> J32(n3, n2);
- gmshMatrix<double> J23(n2, n3);
- gmshVector<double> D3(n3);
- gmshVector<double> R2(n2);
- gmshMatrix<double> J23K33(n2, n3);
+ fullMatrix<double> K33(n3, n3);
+ fullMatrix<double> K22(n2, n2);
+ fullMatrix<double> J32(n3, n2);
+ fullMatrix<double> J23(n2, n3);
+ fullVector<double> D3(n3);
+ fullVector<double> R2(n2);
+ fullMatrix<double> J23K33(n2, n3);
K33.set_all(0.0);
El.elementMatrix(e, K33);
computeMetricVector(gf, e, J32, J23, D3);
@@ -890,7 +890,7 @@ struct swap_triangles_p2
vv.push_back(_test);vv.push_back(n23);vv.push_back(n24);
MTriangleN t4_test (n4,n3,n2,vv,2,t2->getNum(),t2->getPartition());
p2data data(gf,&t3_test,&t4_test,_test,0);
- gmshVector<double> pp(2);
+ fullVector<double> pp(2);
pp(0) = p34_linear.x();
pp(1) = p34_linear.y();
minimize_grad_fd (_function_p2tB, pp, &data);
@@ -1028,7 +1028,7 @@ static int optimalLocationP2_(GFace *gf,
SPoint2 p34_linear = (p1+p2)*.5;
SPoint2 dirt = p4-p3;
- gmshVector<double> pp(2);
+ fullVector<double> pp(2);
pp(0) = p12.x();
pp(1) = p12.y();
p2data data(gf,t1,t2,n12,s);
@@ -1085,7 +1085,7 @@ int optimalLocationPN_ (GFace *gf, const MEdge &me, MTriangle *t1, MTriangle *t2
toOptimize.push_back(t2->getVertex(i));
}
- gmshVector<double> pp(2*toOptimize.size());
+ fullVector<double> pp(2*toOptimize.size());
for (unsigned int i=0;i<toOptimize.size();i++){
SPoint2 pt;
reparamMeshVertexOnFace(toOptimize[i],gf,pt);
diff --git a/Mesh/gmshSmoothHighOrder.h b/Mesh/gmshSmoothHighOrder.h
index 184bd82d0939fe5dab22f455b92724ddff08d6a2..cd536e6219194ec2635031fe1977f339ac05adcb 100644
--- a/Mesh/gmshSmoothHighOrder.h
+++ b/Mesh/gmshSmoothHighOrder.h
@@ -9,7 +9,7 @@
#include <map>
#include <vector>
#include "SVector3.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include "dofManager.h"
#include "elasticityTerm.h"
@@ -56,9 +56,9 @@ public:
faceContainer &faceVertices);
void computeMetricVector(GFace *gf,
MElement *e,
- gmshMatrix<double> &J,
- gmshMatrix<double> &JT,
- gmshVector<double> &D);
+ fullMatrix<double> &J,
+ fullMatrix<double> &JT,
+ fullVector<double> &D);
void moveToStraightSidedLocation(MVertex *) const;
void moveToTargetLocation(MVertex *) const;
void moveToStraightSidedLocation(MElement *) const;
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index aaaee3c204337427aecac2dc950ac49229c959d3..2373738b48d9a2e6fbb145d77ebf455d87300746 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -5,7 +5,7 @@
#include <stdlib.h>
#include "GmshMessage.h"
-#include "GmshPredicates.h"
+#include "robustPredicates.h"
#include "meshGFace.h"
#include "meshGFaceOptimize.h"
#include "BackgroundMesh.h"
@@ -319,8 +319,8 @@ bool edgeSwapTestDelaunay(BDS_Edge *e, GFace *gf)
double op2x[3] = {op[1]->X, op[1]->Y, op[1]->Z};
double fourth[3];
fourthPoint(p1x, p2x, op1x, fourth);
- double result = gmsh::insphere(p1x, p2x, op1x, fourth, op2x) *
- gmsh::orient3d(p1x, p2x, op1x, fourth);
+ double result = robustPredicates::insphere(p1x, p2x, op1x, fourth, op2x) *
+ robustPredicates::orient3d(p1x, p2x, op1x, fourth);
return result > 0.;
}
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index c45149f2150e1dd462a1dc0c70a1678531e6f7f8..20fa4ce22645c2407d3ff0bba8ff506cff8f82d6 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -7,7 +7,7 @@
#include <map>
#include <algorithm>
#include "GmshMessage.h"
-#include "GmshPredicates.h"
+#include "robustPredicates.h"
#include "BackgroundMesh.h"
#include "meshGFaceDelaunayInsertion.h"
#include "meshGFaceOptimize.h"
@@ -89,7 +89,7 @@ void circumCenterMetricXYZ(double *p1, double *p2, double *p3, SMetric3 & metric
double p3P[2]; prosca(v2, vx, &p3P[0]); prosca(v2, vy, &p3P[1]);
double resP[2];
- gmshMatrix<double> T(3,3);
+ fullMatrix<double> T(3,3);
for (int i = 0; i < 3; i++)T(0,i) = vx[i];
for (int i = 0; i < 3; i++)T(1,i) = vy[i];
for (int i = 0; i < 3; i++)T(2,i) = vz[i];
@@ -245,8 +245,8 @@ int MTri3::inCircumCircle(const double *p) const
double fourth[3];
fourthPoint(pa, pb, pc, fourth);
- double result = gmsh::insphere(pa, pb, pc, fourth, (double*)p) *
- gmsh::orient3d(pa, pb, pc,fourth);
+ double result = robustPredicates::insphere(pa, pb, pc, fourth, (double*)p) *
+ robustPredicates::orient3d(pa, pb, pc,fourth);
return (result > 0) ? 1 : 0;
}
@@ -261,8 +261,8 @@ int inCircumCircle(MTriangle *base,
double pc[2] = {Us[base->getVertex(2)->getIndex()],
Vs[base->getVertex(2)->getIndex()]};
- double result =
- gmsh::incircle(pa, pb, pc, (double*)param) * gmsh::orient2d(pa, pb, pc);
+ double result = robustPredicates::incircle(pa, pb, pc, (double*)param) *
+ robustPredicates::orient2d(pa, pb, pc);
return (result > 0) ? 1 : 0;
}
@@ -761,8 +761,8 @@ void testTensor()
t(1,0) = .2;
t(1,1) = 2;
t(2,2) = 3;
- gmshMatrix<double> m(3,3);
- gmshVector<double> v(3);
+ fullMatrix<double> m(3,3);
+ fullVector<double> v(3);
t.eig(m,v);
printf("%12.5E %12.5E %12.5E \n",v(0),v(1),v(2));
printf("%12.5E %12.5E %12.5E \n",m(0,0),m(1,0),m(2,0));
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 96fd06270a47bcbd87cdeb74c9a588d5d11ff382..ade5b52d4842e9b3e53f622cb5e3fb499c1cf6b0 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -7,7 +7,7 @@
#include <map>
#include <algorithm>
#include "GmshMessage.h"
-#include "GmshPredicates.h"
+#include "robustPredicates.h"
#include "OS.h"
#include "BackgroundMesh.h"
#include "meshGRegion.h"
@@ -32,8 +32,8 @@ int MTet4::inCircumSphere(const double *p) const
double pd[3] = {base->getVertex(3)->x(),
base->getVertex(3)->y(),
base->getVertex(3)->z()};
- double result =
- gmsh::insphere(pa, pb, pc, pd, (double*)p) * gmsh::orient3d(pa, pb, pc, pd);
+ double result = robustPredicates::insphere(pa, pb, pc, pd, (double*)p) *
+ robustPredicates::orient3d(pa, pb, pc, pd);
return (result > 0) ? 1 : 0;
}
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index a7eb3c329757072d932a461ea3a52840777d2262..12d2582611696a170c6cb56ab4ddc6fd13cdded9 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -9,7 +9,7 @@
#include "MTriangle.h"
#include "MTetrahedron.h"
#include "Numeric.h"
-#include "FunctionSpace.h"
+#include "functionSpace.h"
#include "GmshMessage.h"
double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3,
@@ -234,7 +234,7 @@ double qmDistorsionOfMapping (MTriangle *e)
const double di = mesh_functional_distorsion (e, u, v);
dmin = (i == 0)? di : std::min(dmin, di);
}
- const gmshMatrix<double>& points = e->getFunctionSpace()->points;
+ const fullMatrix<double>& points = e->getFunctionSpace()->points;
for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
const double u = points(i, 0);
@@ -281,7 +281,7 @@ double qmDistorsionOfMapping(MTetrahedron *e)
dmin = (i == 0) ? di : std::min(dmin, di);
}
- const gmshMatrix<double>& points = e->getFunctionSpace()->points;
+ const fullMatrix<double>& points = e->getFunctionSpace()->points;
for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
const double u = points(i, 0);
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 9824c8e24265151195a9d3b774c8d379f9ccdabc..86deeda39e224be0af7bd9baa20d0914f5a394a0 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -5,26 +5,26 @@
set(SRC
Numeric.cpp
- GmshMatrix.cpp
- gmshLaplace.cpp
- gmshCrossConf.cpp
- gmshDistance.cpp
- gmshConvexCombination.cpp
- gmshHelmholtz.cpp
- gmshElasticity.cpp
- gmshProjection.cpp
- gmshLinearSystemCSR.cpp
+ fullMatrix.cpp
EigSolve.cpp
- FunctionSpace.cpp
- GaussQuadratureLin.cpp
- GaussQuadratureTri.cpp
- GaussQuadratureQuad.cpp
- GaussQuadratureTet.cpp
- GaussQuadratureHex.cpp
- GaussLegendreSimplex.cpp
- GmshPredicates.cpp
+ functionSpace.cpp
+ GaussQuadratureLin.cpp
+ GaussQuadratureTri.cpp
+ GaussQuadratureQuad.cpp
+ GaussQuadratureTet.cpp
+ GaussQuadratureHex.cpp
+ GaussLegendreSimplex.cpp
+ robustPredicates.cpp
mathEvaluator.cpp
Iso.cpp
+ gmshLaplace.cpp
+ gmshCrossConf.cpp
+ gmshDistance.cpp
+ gmshConvexCombination.cpp
+ gmshHelmholtz.cpp
+ gmshElasticity.cpp
+ gmshProjection.cpp
+ gmshLinearSystemCSR.cpp
)
append_gmsh_src(Numeric "${SRC}")
diff --git a/Numeric/FunctionSpace.cpp b/Numeric/FunctionSpace.cpp
deleted file mode 100644
index 59b1b0d63a3ef87c3f4114aa2367089fac0cec21..0000000000000000000000000000000000000000
--- a/Numeric/FunctionSpace.cpp
+++ /dev/null
@@ -1,672 +0,0 @@
-// 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>.
-//
-// Contributor(s):
-// Koen Hillewaert
-//
-
-#include "FunctionSpace.h"
-#include "GmshDefines.h"
-#include "GmshMessage.h"
-
-static gmshMatrix<double> generate1DMonomials(int order)
-{
- gmshMatrix<double> monomials(order + 1, 1);
- for (int i = 0; i < order + 1; i++) monomials(i, 0) = i;
- return monomials;
-}
-
-static gmshMatrix<double> generate1DPoints(int order)
-{
- gmshMatrix<double> line(order + 1, 1);
- line(0, 0) = -1.;
- line(1, 0) = 1.;
- double dd = 2. / order;
- for (int i = 2; i < order + 1; i++) line(i, 0) = -1. + dd * (i - 1);
- return line;
-}
-
-static gmshMatrix<double> generatePascalTriangle(int order)
-{
- gmshMatrix<double> monomials((order + 1) * (order + 2) / 2, 2);
- int index = 0;
- for (int i = 0; i <= order; i++) {
- for (int j = 0; j <= i; j++) {
- monomials(index, 0) = i - j;
- monomials(index, 1) = j;
- index++;
- }
- }
- return monomials;
-}
-
-// generate the exterior hull of the Pascal triangle for Serendipity element
-
-static gmshMatrix<double> generatePascalSerendipityTriangle(int order)
-{
- gmshMatrix<double> monomials(3 * order, 2);
-
- monomials(0, 0) = 0;
- monomials(0, 1) = 0;
-
- int index = 1;
- for (int i = 1; i <= order; i++) {
- if (i == order) {
- for (int j = 0; j <= i; j++) {
- monomials(index, 0) = i - j;
- monomials(index, 1) = j;
- index++;
- }
- }
- else {
- monomials(index, 0) = i;
- monomials(index, 1) = 0;
- index++;
- monomials(index, 0) = 0;
- monomials(index, 1) = i;
- index++;
- }
- }
- return monomials;
-}
-
-// generate the monomials subspace of all monomials of order exactly == p
-
-static gmshMatrix<double> generateMonomialSubspace(int dim, int p)
-{
- gmshMatrix<double> monomials;
-
- switch (dim) {
- case 1:
- monomials = gmshMatrix<double>(1, 1);
- monomials(0, 0) = p;
- break;
- case 2:
- monomials = gmshMatrix<double>(p + 1, 2);
- for (int k = 0; k <= p; k++) {
- monomials(k, 0) = p - k;
- monomials(k, 1) = k;
- }
- break;
- case 3:
- monomials = gmshMatrix<double>((p + 1) * (p + 2) / 2, 3);
- int index = 0;
- for (int i = 0; i <= p; i++) {
- for (int k = 0; k <= p - i; k++) {
- monomials(index, 0) = p - i - k;
- monomials(index, 1) = k;
- monomials(index, 2) = i;
- index++;
- }
- }
- break;
- }
- return monomials;
-}
-
-// generate external hull of the Pascal tetrahedron
-
-static gmshMatrix<double> generatePascalSerendipityTetrahedron(int order)
-{
- int nbMonomials = 4 + 6 * std::max(0, order - 1) +
- 4 * std::max(0, (order - 2) * (order - 1) / 2);
- gmshMatrix<double> monomials(nbMonomials, 3);
-
- monomials.set_all(0);
- int index = 1;
- for (int p = 1; p < order; p++) {
- for (int i = 0; i < 3; i++) {
- int j = (i + 1) % 3;
- int k = (i + 2) % 3;
- for (int ii = 0; ii < p; ii++) {
- monomials(index, i) = p - ii;
- monomials(index, j) = ii;
- monomials(index, k) = 0;
- index++;
- }
- }
- }
- gmshMatrix<double> monomialsMaxOrder = generateMonomialSubspace(3, order);
- int nbMaxOrder = monomialsMaxOrder.size1();
- monomials.copy(monomialsMaxOrder, 0, nbMaxOrder, 0, 3, index, 0);
- return monomials;
-}
-
-// generate Pascal tetrahedron
-
-static gmshMatrix<double> generatePascalTetrahedron(int order)
-{
- int nbMonomials = (order + 1) * (order + 2) * (order + 3) / 6;
-
- gmshMatrix<double> monomials(nbMonomials, 3);
-
- int index = 0;
- for (int p = 0; p <= order; p++) {
- gmshMatrix<double> monOrder = generateMonomialSubspace(3, p);
- int nb = monOrder.size1();
- monomials.copy(monOrder, 0, nb, 0, 3, index, 0);
- index += nb;
- }
-
- return monomials;
-}
-
-static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
-//static int nbdoftriangleserendip(int order) { return 3 * order; }
-
-//KH : caveat : node coordinates are not yet coherent with node numbering associated
-// to numbering of principal vertices of face !!!!
-
-// uv surface - orientation v0-v2-v1
-static void nodepositionface0(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= order; w[1] = 0.;
- u[2]= order; v[2]= 0.; w[2] = 0.;
-
- // edges
- for (int k = 0; k < (order - 1); k++){
- u[3 + k] = 0.;
- v[3 + k] = k + 1;
- w[3 + k] = 0.;
-
- u[3 + order - 1 + k] = k + 1;
- v[3 + order - 1 + k] = order - 1 - k ;
- w[3 + order - 1 + k] = 0.;
-
- u[3 + 2 * (order - 1) + k] = order - 1 - k;
- v[3 + 2 * (order - 1) + k] = 0.;
- w[3 + 2 * (order - 1) + k] = 0.;
- }
-
- if (order > 2){
- int nbdoftemp = nbdoftriangle(order - 3);
- nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
- &w[3 + 3* (order - 1)]);
- for (int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3);
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-// uw surface - orientation v0-v1-v3
-static void nodepositionface1(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
- u[0] = 0.; v[0]= 0.; w[0] = 0.;
- u[1] = order; v[1]= 0.; w[1] = 0.;
- u[2] = 0.; v[2]= 0.; w[2] = order;
- // edges
- for (int k = 0; k < (order - 1); k++){
- u[3 + k] = k + 1;
- v[3 + k] = 0.;
- w[3 + k] = 0.;
-
- u[3 + order - 1 + k] = order - 1 - k;
- v[3 + order - 1 + k] = 0.;
- w[3 + order - 1+ k ] = k + 1;
-
- u[3 + 2 * (order - 1) + k] = 0. ;
- v[3 + 2 * (order - 1) + k] = 0.;
- w[3 + 2 * (order - 1) + k] = order - 1 - k;
- }
- if (order > 2){
- int nbdoftemp = nbdoftriangle(order - 3);
- nodepositionface1(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order -1 )],
- &w[3 + 3 * (order - 1)]);
- for (int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-// vw surface - orientation v0-v3-v2
-static void nodepositionface2(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
-
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= 0.; w[1] = order;
- u[2]= 0.; v[2]= order; w[2] = 0.;
- // edges
- for (int k = 0; k < (order - 1); k++){
-
- u[3 + k] = 0.;
- v[3 + k] = 0.;
- w[3 + k] = k + 1;
-
- u[3 + order - 1 + k] = 0.;
- v[3 + order - 1 + k] = k + 1;
- w[3 + order - 1 + k] = order - 1 - k;
-
- u[3 + 2 * (order - 1) + k] = 0.;
- v[3 + 2 * (order - 1) + k] = order - 1 - k;
- w[3 + 2 * (order - 1) + k] = 0.;
- }
- if (order > 2){
- int nbdoftemp = nbdoftriangle(order - 3);
- nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
- &w[3 + 3 * (order - 1)]);
- for (int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-// uvw surface - orientation v3-v1-v2
-static void nodepositionface3(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
- u[0]= 0.; v[0]= 0.; w[0] = order;
- u[1]= order; v[1]= 0.; w[1] = 0.;
- u[2]= 0.; v[2]= order; w[2] = 0.;
- // edges
- for (int k = 0; k < (order - 1); k++){
-
- u[3 + k] = k + 1;
- v[3 + k] = 0.;
- w[3 + k] = order - 1 - k;
-
- u[3 + order - 1 + k] = order - 1 - k;
- v[3 + order - 1 + k] = k + 1;
- w[3 + order - 1 + k] = 0.;
-
- u[3 + 2 * (order - 1) + k] = 0.;
- v[3 + 2 * (order - 1) + k] = order - 1 - k;
- w[3 + 2 * (order - 1) + k] = k + 1;
- }
- if (order > 2){
- int nbdoftemp = nbdoftriangle(order - 3);
- nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
- &w[3 + 3 * (order - 1)]);
- for (int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-static gmshMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip)
-{
- int nbPoints =
- (serendip ?
- 4 + 6 * std::max(0, order - 1) + 4 * std::max(0, (order - 2) * (order - 1) / 2) :
- (order + 1) * (order + 2) * (order + 3) / 6);
-
- gmshMatrix<double> point(nbPoints, 3);
-
- double overOrder = (order == 0 ? 1. : 1. / order);
-
- point(0, 0) = 0.;
- point(0, 1) = 0.;
- point(0, 2) = 0.;
-
- if (order > 0) {
- point(1, 0) = order;
- point(1, 1) = 0;
- point(1, 2) = 0;
-
- point(2, 0) = 0.;
- point(2, 1) = order;
- point(2, 2) = 0.;
-
- point(3, 0) = 0.;
- point(3, 1) = 0.;
- point(3, 2) = order;
-
- // edges e5 and e6 switched in original version, opposite direction
- // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
-
- if (order > 1) {
- for (int k = 0; k < (order - 1); k++) {
- point(4 + k, 0) = k + 1;
- point(4 + order - 1 + k, 0) = order - 1 - k;
- point(4 + 2 * (order - 1) + k, 0) = 0.;
- point(4 + 3 * (order - 1) + k, 0) = 0.;
- // point(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
- // point(4 + 5 * (order - 1) + k, 0) = 0.;
- point(4 + 4 * (order - 1) + k, 0) = 0.;
- point(4 + 5 * (order - 1) + k, 0) = k+1;
-
- point(4 + k, 1) = 0.;
- point(4 + order - 1 + k, 1) = k + 1;
- point(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
- point(4 + 3 * (order - 1) + k, 1) = 0.;
- // point(4 + 4 * (order - 1) + k, 1) = 0.;
- // point(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
- point(4 + 4 * (order - 1) + k, 1) = k+1;
- point(4 + 5 * (order - 1) + k, 1) = 0.;
-
- point(4 + k, 2) = 0.;
- point(4 + order - 1 + k, 2) = 0.;
- point(4 + 2 * (order - 1) + k, 2) = 0.;
- point(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
- point(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
- point(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
- }
-
- if (order > 2) {
- int ns = 4 + 6 * (order - 1);
- int nbdofface = nbdoftriangle(order - 3);
-
- double *u = new double[nbdofface];
- double *v = new double[nbdofface];
- double *w = new double[nbdofface];
-
- nodepositionface0(order - 3, u, v, w);
-
- // u-v plane
-
- for (int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3) + 1.;
- point(ns + i, 1) = v[i] * (order - 3) + 1.;
- point(ns + i, 2) = w[i] * (order - 3);
- }
-
- ns = ns + nbdofface;
-
- // u-w plane
-
- nodepositionface1(order - 3, u, v, w);
-
- for (int i=0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3) + 1.;
- point(ns + i, 1) = v[i] * (order - 3) ;
- point(ns + i, 2) = w[i] * (order - 3) + 1.;
- }
-
- // v-w plane
-
- ns = ns + nbdofface;
-
- nodepositionface2(order - 3, u, v, w);
-
- for (int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3);
- point(ns + i, 1) = v[i] * (order - 3) + 1.;
- point(ns + i, 2) = w[i] * (order - 3) + 1.;
- }
-
- // u-v-w plane
-
- ns = ns + nbdofface;
-
- nodepositionface3(order - 3, u, v, w);
-
- for (int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3) + 1.;
- point(ns + i, 1) = v[i] * (order - 3) + 1.;
- point(ns + i, 2) = w[i] * (order - 3) + 1.;
- }
-
- ns = ns + nbdofface;
-
- delete [] u;
- delete [] v;
- delete [] w;
-
- if (!serendip && order > 3) {
-
- gmshMatrix<double> interior = gmshGeneratePointsTetrahedron(order - 4, false);
- for (int k = 0; k < interior.size1(); k++) {
- point(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
- point(ns + k, 1) = 1. + interior(k, 1) * (order - 4);
- point(ns + k, 2) = 1. + interior(k, 2) * (order - 4);
- }
- }
- }
- }
- }
-
- point.scale(overOrder);
- return point;
-}
-
-static gmshMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
-{
- int nbPoints = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
- gmshMatrix<double> point(nbPoints, 2);
-
- point(0, 0) = 0;
- point(0, 1) = 0;
-
- double dd = 1. / order;
-
- if (order > 0) {
- point(1, 0) = 1;
- point(1, 1) = 0;
- point(2, 0) = 0;
- point(2, 1) = 1;
-
- int index = 3;
-
- if (order > 1) {
-
- double ksi = 0;
- double eta = 0;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi += dd;
- point(index, 0) = ksi;
- point(index, 1) = eta;
- }
-
- ksi = 1.;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi -= dd;
- eta += dd;
- point(index, 0) = ksi;
- point(index, 1) = eta;
- }
-
- eta = 1.;
- ksi = 0.;
-
- for (int i = 0; i < order - 1; i++, index++) {
- eta -= dd;
- point(index, 0) = ksi;
- point(index, 1) = eta;
- }
-
- if (order > 2 && !serendip) {
- gmshMatrix<double> inner = gmshGeneratePointsTriangle(order - 3, serendip);
- inner.scale(1. - 3. * dd);
- inner.add(dd);
- point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
- return point;
-}
-
-static gmshMatrix<double> generateLagrangeMonomialCoefficients
- (const gmshMatrix<double>& monomial, const gmshMatrix<double>& point)
-{
- if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
- Msg::Fatal("Wrong sizes for Lagrange coefficients generation");
- return gmshMatrix<double>(1, 1);
- }
-
- int ndofs = monomial.size1();
- int dim = monomial.size2();
-
- gmshMatrix<double> Vandermonde(ndofs, ndofs);
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- double dd = 1.;
- for (int k = 0; k < dim; k++) dd *= pow(point(j, k), monomial(i, k));
- Vandermonde(i, j) = dd;
- }
- }
-
- gmshMatrix<double> coefficient(ndofs, ndofs);
- Vandermonde.invert(coefficient);
- return coefficient;
-}
-
-std::map<int, gmshFunctionSpace> gmshFunctionSpaces::fs;
-
-const gmshFunctionSpace &gmshFunctionSpaces::find(int tag)
-{
- std::map<int, gmshFunctionSpace>::const_iterator it = fs.find(tag);
- if (it != fs.end()) return it->second;
-
- gmshFunctionSpace F;
-
- switch (tag){
- case MSH_PNT:
- F.monomials = generate1DMonomials(0);
- F.points = generate1DPoints(0);
- break;
- case MSH_LIN_2 :
- F.monomials = generate1DMonomials(1);
- F.points = generate1DPoints(1);
- break;
- case MSH_LIN_3 :
- F.monomials = generate1DMonomials(2);
- F.points = generate1DPoints(2);
- break;
- case MSH_LIN_4:
- F.monomials = generate1DMonomials(3);
- F.points = generate1DPoints(3);
- break;
- case MSH_LIN_5:
- F.monomials = generate1DMonomials(4);
- F.points = generate1DPoints(4);
- break;
- case MSH_LIN_6:
- F.monomials = generate1DMonomials(5);
- F.points = generate1DPoints(5);
- break;
- case MSH_TRI_3 :
- F.monomials = generatePascalTriangle(1);
- F.points = gmshGeneratePointsTriangle(1, false);
- break;
- case MSH_TRI_6 :
- F.monomials = generatePascalTriangle(2);
- F.points = gmshGeneratePointsTriangle(2, false);
- break;
- case MSH_TRI_9 :
- F.monomials = generatePascalSerendipityTriangle(3);
- F.points = gmshGeneratePointsTriangle(3, true);
- break;
- case MSH_TRI_10 :
- F.monomials = generatePascalTriangle(3);
- F.points = gmshGeneratePointsTriangle(3, false);
- break;
- case MSH_TRI_12 :
- F.monomials = generatePascalSerendipityTriangle(4);
- F.points = gmshGeneratePointsTriangle(4, true);
- break;
- case MSH_TRI_15 :
- F.monomials = generatePascalTriangle(4);
- F.points = gmshGeneratePointsTriangle(4, false);
- break;
- case MSH_TRI_15I :
- F.monomials = generatePascalSerendipityTriangle(5);
- F.points = gmshGeneratePointsTriangle(5, true);
- break;
- case MSH_TRI_21 :
- F.monomials = generatePascalTriangle(5);
- F.points = gmshGeneratePointsTriangle(5, false);
- break;
- case MSH_TET_4 :
- F.monomials = generatePascalTetrahedron(1);
- F.points = gmshGeneratePointsTetrahedron(1, false);
- break;
- case MSH_TET_10 :
- F.monomials = generatePascalTetrahedron(2);
- F.points = gmshGeneratePointsTetrahedron(2, false);
- break;
- case MSH_TET_20 :
- F.monomials = generatePascalTetrahedron(3);
- F.points = gmshGeneratePointsTetrahedron(3, false);
- break;
- case MSH_TET_35 :
- F.monomials = generatePascalTetrahedron(4);
- F.points = gmshGeneratePointsTetrahedron(4, false);
- break;
- case MSH_TET_34 :
- F.monomials = generatePascalSerendipityTetrahedron(4);
- F.points = gmshGeneratePointsTetrahedron(4, true);
- break;
- case MSH_TET_52 :
- F.monomials = generatePascalSerendipityTetrahedron(5);
- F.points = gmshGeneratePointsTetrahedron(5, true);
- break;
- case MSH_TET_56 :
- F.monomials = generatePascalTetrahedron(5);
- F.points = gmshGeneratePointsTetrahedron(5, false);
- break;
- default :
- Msg::Error("Unknown function space %d: reverting to TET_4", tag);
- F.monomials = generatePascalTetrahedron(1);
- F.points = gmshGeneratePointsTetrahedron(1, false);
- break;
- }
- F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
- fs.insert(std::make_pair(tag, F));
- return fs[tag];
-}
-
-std::map<std::pair<int, int>, gmshMatrix<double> > gmshFunctionSpaces::injector;
-
-const gmshMatrix<double> &gmshFunctionSpaces::findInjector(int tag1, int tag2)
-{
- std::pair<int,int> key(tag1,tag2);
- std::map<std::pair<int, int>, gmshMatrix<double> >::const_iterator it = injector.find(key);
- if (it != injector.end()) return it->second;
-
- const gmshFunctionSpace& fs1 = find(tag1);
- const gmshFunctionSpace& fs2 = find(tag2);
-
- gmshMatrix<double> inj(fs1.points.size1(),fs2.points.size1());
-
- double sf[256];
-
- for (int i = 0; i < fs1.points.size1(); i++) {
- fs2.f(fs1.points(i, 0), fs1.points(i, 1), fs1.points(i, 2), sf);
- for (int j = 0; j < fs2.points.size1(); j++) inj(i, j) = sf[j];
- }
-
- injector.insert(std::make_pair(key, inj));
- return injector[key];
-}
diff --git a/Numeric/FunctionSpace.h b/Numeric/FunctionSpace.h
deleted file mode 100644
index d701cf2fadc5e95dfb2dc88b59b6b28bc3919e78..0000000000000000000000000000000000000000
--- a/Numeric/FunctionSpace.h
+++ /dev/null
@@ -1,107 +0,0 @@
-// 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 _FUNCTION_SPACE_H_
-#define _FUNCTION_SPACE_H_
-
-#include <math.h>
-#include <map>
-#include "GmshMatrix.h"
-
-struct gmshFunctionSpace
-{
- gmshMatrix<double> points;
- gmshMatrix<double> monomials;
- gmshMatrix<double> coefficients;
- inline void evaluateMonomials(double u, double v, double w, double p[]) const
- {
- for (int j = 0; j < monomials.size1(); j++) {
- p[j] = pow(u, (int)monomials(j, 0));
- if (monomials.size2() > 1) p[j] *= pow(v, (int)monomials(j, 1));
- if (monomials.size2() > 2) p[j] *= pow(w, (int)monomials(j, 2));
- }
- }
- inline void f(double u, double v, double w, double *sf) const
- {
- double p[256];
- evaluateMonomials(u, v, w, p);
- for (int i = 0; i < coefficients.size1(); i++) {
- sf[i] = 0;
- for (int j = 0; j < coefficients.size2(); j++) {
- sf[i] += coefficients(i, j) * p[j];
- }
- }
- }
- inline void df(double u, double v, double w, double grads[][3]) const
- {
- switch (monomials.size2()) {
- case 1:
- for (int i = 0; i < coefficients.size1(); i++){
- grads[i][0] = 0;
- grads[i][1] = 0;
- grads[i][2] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- if ((monomials)(j, 0) > 0)
- grads[i][0] += (coefficients)(i, j) *
- pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0);
- }
- }
- break;
- case 2:
- for (int i = 0; i < coefficients.size1(); i++){
- grads[i][0] = 0;
- grads[i][1] = 0;
- grads[i][2] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- if ((monomials)(j, 0) > 0)
- grads[i][0] += (coefficients)(i, j) *
- pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) *
- pow(v, (monomials)(j, 1));
- if ((monomials)(j, 1) > 0)
- grads[i][1] += (coefficients)(i, j) *
- pow(u, (monomials)(j, 0)) *
- pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1);
- }
- }
- break;
- case 3:
- for (int i = 0; i < coefficients.size1(); i++){
- grads[i][0] = 0;
- grads[i][1] = 0;
- grads[i][2] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- if ((monomials)(j, 0) > 0)
- grads[i][0] += (coefficients)(i, j) *
- pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) *
- pow(v, (monomials)(j, 1)) *
- pow(w, (monomials)(j, 2));
- if ((monomials)(j, 1) > 0)
- grads[i][1] += (coefficients)(i, j) *
- pow(u, (monomials)(j, 0)) *
- pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1) *
- pow(w, (monomials)(j, 2));
- if ((monomials)(j, 2) > 0)
- grads[i][2] += (coefficients)(i, j) *
- pow(u, (monomials)(j, 0)) *
- pow(v, (monomials)(j, 1)) *
- pow(w, (monomials)(j, 2) - 1) * (monomials)(j, 2);
- }
- }
- break;
- }
- }
-};
-
-class gmshFunctionSpaces
-{
- private:
- static std::map<int, gmshFunctionSpace> fs;
- static std::map<std::pair<int, int>, gmshMatrix<double> > injector;
- public :
- static const gmshFunctionSpace &find(int);
- static const gmshMatrix<double> &findInjector(int, int);
-};
-
-#endif
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 8285b1df9850b8a5159157b122cb7aa9980040fe..fe6fcfb72d0087f5e898d0f7f3d4616d381fd13e 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -511,8 +511,8 @@ void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
}
}
- gmshMatrix<double> M(3, 3), V(3, 3);
- gmshVector<double> W(3);
+ fullMatrix<double> M(3, 3), V(3, 3);
+ fullVector<double> W(3);
for(i = 1; i <= n; i++){
for(j = 1; j <= n; j++){
M(i - 1, j - 1) = MM[i - 1][j - 1];
@@ -536,15 +536,15 @@ void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
}
}
-bool newton_fd(void (*func)(gmshVector<double> &, gmshVector<double> &, void *),
- gmshVector<double> &x, void *data, double relax, double tolx)
+bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
+ fullVector<double> &x, void *data, double relax, double tolx)
{
const int MAXIT = 50;
const double EPS = 1.e-4;
const int N = x.size();
- gmshMatrix<double> J(N, N);
- gmshVector<double> f(N), feps(N), dx(N);
+ fullMatrix<double> J(N, N);
+ fullVector<double> f(N), feps(N), dx(N);
for (int iter = 0; iter < MAXIT; iter++){
func(x, f, data);
@@ -588,9 +588,9 @@ f(x+a*d) = f(x) + f'(x) (
*/
-void gmshLineSearch(double (*func)(gmshVector<double> &, void *), void* data,
- gmshVector<double> &x, gmshVector<double> &p,
- gmshVector<double> &g, double &f,
+void gmshLineSearch(double (*func)(fullVector<double> &, void *), void* data,
+ fullVector<double> &x, fullVector<double> &p,
+ fullVector<double> &g, double &f,
double stpmax, int &check)
{
int i;
@@ -599,7 +599,7 @@ void gmshLineSearch(double (*func)(gmshVector<double> &, void *), void* data,
const double ALF = 1.0e-4;
const double TOLX = 1.0e-9;
- gmshVector<double> xold(x);
+ fullVector<double> xold(x);
const double fold = (*func)(xold,data);
check=0;
@@ -659,15 +659,15 @@ void gmshLineSearch(double (*func)(gmshVector<double> &, void *), void* data,
}
}
-double minimize_grad_fd (double (*func)(gmshVector<double> &, void *),
- gmshVector<double> &x, void *data)
+double minimize_grad_fd (double (*func)(fullVector<double> &, void *),
+ fullVector<double> &x, void *data)
{
const int MAXIT = 3;
const double EPS = 1.e-4;
const int N = x.size();
- gmshVector<double> grad(N);
- gmshVector<double> dir(N);
+ fullVector<double> grad(N);
+ fullVector<double> dir(N);
double f,feps,finit;
for (int iter = 0; iter < MAXIT; iter++){
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index fd672af9cd16ea130b2189a6e2383e8fc06c99c0..134be5b0b15f1c36dd576a497ae38dfbc773b85e 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -7,7 +7,7 @@
#define _NUMERIC_H_
#include <math.h>
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#define myhypot(a,b) (sqrt((a)*(a)+(b)*(b)))
#define sign(x) (((x)>=0)?1:-1)
@@ -72,9 +72,9 @@ void gradSimplex(double *x, double *y, double *z, double *v, double *grad);
double ComputeVonMises(double *val);
double ComputeScalarRep(int numComp, double *val);
void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
-bool newton_fd(void (*func)(gmshVector<double> &, gmshVector<double> &, void *),
- gmshVector<double> &x, void *data, double relax=1., double tolx=1.e-6);
-double minimize_grad_fd (double (*func)(gmshVector<double> &, void *),
- gmshVector<double> &x, void *data);
+bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
+ fullVector<double> &x, void *data, double relax=1., double tolx=1.e-6);
+double minimize_grad_fd(double (*func)(fullVector<double> &, void *),
+ fullVector<double> &x, void *data);
#endif
diff --git a/Numeric/gmshConvexCombination.cpp b/Numeric/gmshConvexCombination.cpp
index 814bb2a3189ddcf5bcd1178c2d7495e8b9158150..8f0f2d9315a513c42589c9c97abb3b7a61576656 100644
--- a/Numeric/gmshConvexCombination.cpp
+++ b/Numeric/gmshConvexCombination.cpp
@@ -8,7 +8,7 @@
#include "gmshConvexCombination.h"
#include "Numeric.h"
-void gmshConvexCombinationTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+void gmshConvexCombinationTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
{
m.set_all(0.);
int nbNodes = e->getNumVertices();
diff --git a/Numeric/gmshConvexCombination.h b/Numeric/gmshConvexCombination.h
index a2ae7c321bb1bde1924461bfce3186aaa838da64..bb946e41d3668204b20be2b7361c3fc2881c357b 100644
--- a/Numeric/gmshConvexCombination.h
+++ b/Numeric/gmshConvexCombination.h
@@ -11,7 +11,7 @@
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
class gmshConvexCombinationTerm : public gmshNodalFemTerm<double> {
protected:
@@ -29,7 +29,7 @@ class gmshConvexCombinationTerm : public gmshNodalFemTerm<double> {
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;
+ virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
};
diff --git a/Numeric/gmshCrossConf.cpp b/Numeric/gmshCrossConf.cpp
index 938bafbf0c7410df7b7e20033850e52b0a4b06ec..f81d5d3562152e0b0f33f6f0b121fc0c6db0a894 100644
--- a/Numeric/gmshCrossConf.cpp
+++ b/Numeric/gmshCrossConf.cpp
@@ -8,7 +8,7 @@
#include "gmshCrossConf.h"
#include "Numeric.h"
-void gmshCrossConfTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+void gmshCrossConfTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
diff --git a/Numeric/gmshCrossConf.h b/Numeric/gmshCrossConf.h
index 801b52ba084697e6c4ba19e9899eb6ae0cc22db0..792aac05ec6b566315e2b76d5016e67ef519f88a 100644
--- a/Numeric/gmshCrossConf.h
+++ b/Numeric/gmshCrossConf.h
@@ -13,7 +13,7 @@
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
class gmshCrossConfTerm : public gmshNodalFemTerm<double> {
protected:
@@ -41,7 +41,7 @@ class gmshCrossConfTerm : public gmshNodalFemTerm<double> {
{
_iFieldB = (iFieldB==-1) ? _iField : iFieldB;
}
- virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+ virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
};
class gmshCrossConfTerm2DParametric : gmshCrossConfTerm {
@@ -49,7 +49,7 @@ class gmshCrossConfTerm2DParametric : gmshCrossConfTerm {
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;
+ virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
};
diff --git a/Numeric/gmshDistance.cpp b/Numeric/gmshDistance.cpp
index ec11b25abc954623a8c050f48628ab6db6e8decc..506d777210169c32e320402f8f2959e6f01b1732 100644
--- a/Numeric/gmshDistance.cpp
+++ b/Numeric/gmshDistance.cpp
@@ -8,7 +8,7 @@
#include "gmshDistance.h"
#include "Numeric.h"
-void gmshDistanceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+void gmshDistanceTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
@@ -58,7 +58,7 @@ void gmshDistanceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
}
-void gmshDistanceTerm::elementVector(MElement *e, gmshVector<double> &m) const{
+void gmshDistanceTerm::elementVector(MElement *e, fullVector<double> &m) const{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
diff --git a/Numeric/gmshDistance.h b/Numeric/gmshDistance.h
index 588cdaec207a048144ad3ac738ede4b1631f6517..57952ea7d4afc20bfc3d1a541e1fd1dfd6c2e1cf 100644
--- a/Numeric/gmshDistance.h
+++ b/Numeric/gmshDistance.h
@@ -13,7 +13,7 @@
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
class gmshDistanceTerm : public gmshNodalFemTerm<double> {
protected:
@@ -41,8 +41,8 @@ class gmshDistanceTerm : public gmshNodalFemTerm<double> {
{
_iFieldB = (iFieldB==-1) ? _iField : iFieldB;
}
- void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
- void elementVector(MElement *e, gmshVector<double> &m) const;
+ void elementMatrix(MElement *e, fullMatrix<double> &m) const;
+ void elementVector(MElement *e, fullVector<double> &m) const;
};
diff --git a/Numeric/gmshElasticity.cpp b/Numeric/gmshElasticity.cpp
index f9e1990283d925d2ae93548f69d4218b321deb6a..1acd57bda9d5db0bf3abf1cb9b8eefafb118d61c 100644
--- a/Numeric/gmshElasticity.cpp
+++ b/Numeric/gmshElasticity.cpp
@@ -14,7 +14,7 @@ Non linear version :
*/
-void gmshElasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+void gmshElasticityTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
@@ -39,10 +39,10 @@ void gmshElasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
{ 0, 0, 0, 0, C44, 0},
{ 0, 0, 0, 0, 0, C44} };
- gmshMatrix<double> H(6, 6);
- gmshMatrix<double> B(6, 3 * nbNodes);
- gmshMatrix<double> BTH(3 * nbNodes, 6);
- gmshMatrix<double> BT(3 * nbNodes, 6);
+ fullMatrix<double> H(6, 6);
+ fullMatrix<double> B(6, 3 * nbNodes);
+ fullMatrix<double> BTH(3 * nbNodes, 6);
+ fullMatrix<double> BT(3 * nbNodes, 6);
for (int i = 0; i < 6; i++)
for (int j = 0; j < 6; j++)
H(i, j) = C[i][j];
@@ -83,7 +83,7 @@ void gmshElasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
}
}
-void gmshElasticityTerm::elementVector(MElement *e, gmshVector<double> &m) const{
+void gmshElasticityTerm::elementVector(MElement *e, fullVector<double> &m) const{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
int npts;
diff --git a/Numeric/gmshElasticity.h b/Numeric/gmshElasticity.h
index 7c291ac53fa1d5ec0528d27c5422954aff778c20..bb736aefc1b36180d9b3b10aa59332fa50e46914 100644
--- a/Numeric/gmshElasticity.h
+++ b/Numeric/gmshElasticity.h
@@ -12,7 +12,7 @@
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
class gmshElasticityTerm : public gmshNodalFemTerm<double> {
protected:
@@ -33,8 +33,8 @@ class gmshElasticityTerm : public gmshNodalFemTerm<double> {
gmshElasticityTerm(GModel *gm, double E, double nu, int iField = 1) :
gmshNodalFemTerm<double>(gm), _E(E), _nu(nu), _iField(iField){}
void setVector(const SVector3 &f) {_f = f;}
- void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
- void elementVector(MElement *e, gmshVector<double> &m) const;
+ void elementMatrix(MElement *e, fullMatrix<double> &m) const;
+ void elementVector(MElement *e, fullVector<double> &m) const;
};
#endif
diff --git a/Numeric/gmshHelmholtz.cpp b/Numeric/gmshHelmholtz.cpp
index 661f40ac7c1d3810c22f7d4f726e56ccde01ff02..1cfb91b793979f3bd03247cdad62387edd74dc5b 100644
--- a/Numeric/gmshHelmholtz.cpp
+++ b/Numeric/gmshHelmholtz.cpp
@@ -9,7 +9,7 @@
#include "Numeric.h"
void gmshHelmholtzTerm::elementMatrix(MElement *e,
- gmshMatrix<std::complex<double> > &m) const
+ fullMatrix<std::complex<double> > &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
diff --git a/Numeric/gmshHelmholtz.h b/Numeric/gmshHelmholtz.h
index 5e08e804d8ac5b173081e63814df32b339c1dc84..800f5416a784d8e7dafe8ed264e35a68cd622daf 100644
--- a/Numeric/gmshHelmholtz.h
+++ b/Numeric/gmshHelmholtz.h
@@ -14,7 +14,7 @@
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
class gmshHelmholtzTerm : public gmshNodalFemTerm<std::complex<double> > {
private:
@@ -33,7 +33,7 @@ class gmshHelmholtzTerm : public gmshNodalFemTerm<std::complex<double> > {
gmshHelmholtzTerm(GModel *gm, gmshFunction<std::complex<double> > *waveNumber,
int iField = 0)
: gmshNodalFemTerm<std::complex<double> >(gm), _waveNumber(waveNumber), _iField(iField){}
- virtual void elementMatrix(MElement *e, gmshMatrix<std::complex<double> > &m) const;
+ virtual void elementMatrix(MElement *e, fullMatrix<std::complex<double> > &m) const;
};
#endif
diff --git a/Numeric/gmshLaplace.cpp b/Numeric/gmshLaplace.cpp
index 83fc8daa913be8b31ee7a695c1429bea79b8ffee..122b3012115e7fd022b6d7d82b5818808c65b4fa 100644
--- a/Numeric/gmshLaplace.cpp
+++ b/Numeric/gmshLaplace.cpp
@@ -8,7 +8,7 @@
#include "gmshLaplace.h"
#include "Numeric.h"
-void gmshLaplaceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+void gmshLaplaceTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
index ae63420ddf7cf82c5c17fbc197f0799b4a0a9026..45c84f429fecdae55f7dcc60391808593e3ce5bb 100644
--- a/Numeric/gmshLaplace.h
+++ b/Numeric/gmshLaplace.h
@@ -13,7 +13,7 @@
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
protected:
@@ -41,7 +41,7 @@ class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
{
_iFieldB = (iFieldB==-1) ? _iField : iFieldB;
}
- virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+ virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
};
diff --git a/Numeric/gmshLinearSystemFull.h b/Numeric/gmshLinearSystemFull.h
index a3c8b129efa72663e8badf93621d99cc2ead45ae..55059b7a2eda81113de9a30dd8385744cabe2d0e 100644
--- a/Numeric/gmshLinearSystemFull.h
+++ b/Numeric/gmshLinearSystemFull.h
@@ -12,24 +12,24 @@
#include "GmshMessage.h"
#include "gmshLinearSystem.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include <stdlib.h>
#include <set>
template <class scalar>
class gmshLinearSystemFull : public gmshLinearSystem<scalar> {
private:
- gmshMatrix<scalar> *_a;
- gmshVector<scalar> *_b, *_x;
+ fullMatrix<scalar> *_a;
+ fullVector<scalar> *_b, *_x;
public :
gmshLinearSystemFull() : _a(0), _b(0), _x(0){}
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int _nbRows)
{
clear();
- _a = new gmshMatrix<scalar>(_nbRows, _nbRows);
- _b = new gmshVector<scalar>(_nbRows);
- _x = new gmshVector<scalar>(_nbRows);
+ _a = new fullMatrix<scalar>(_nbRows, _nbRows);
+ _b = new fullVector<scalar>(_nbRows);
+ _x = new fullVector<scalar>(_nbRows);
}
virtual ~gmshLinearSystemFull()
diff --git a/Numeric/gmshProjection.cpp b/Numeric/gmshProjection.cpp
index 5875c8e668df017e1c02be08e9f1868d0685cb2c..47fb2e159f671ef43f5e185a822f9d87632e16cd 100644
--- a/Numeric/gmshProjection.cpp
+++ b/Numeric/gmshProjection.cpp
@@ -8,7 +8,7 @@
#include "gmshProjection.h"
#include "Numeric.h"
-void gmshProjectionTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+void gmshProjectionTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
diff --git a/Numeric/gmshProjection.h b/Numeric/gmshProjection.h
index d47435ae912a0e8eb7491ee4b1e007cfb00239a6..71287a8b6e982e51d23e5b237e92368dac1f2733 100644
--- a/Numeric/gmshProjection.h
+++ b/Numeric/gmshProjection.h
@@ -13,7 +13,7 @@
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
class gmshProjectionTerm : public gmshNodalFemTerm<double> {
private:
@@ -30,7 +30,7 @@ class gmshProjectionTerm : public gmshNodalFemTerm<double> {
public:
gmshProjectionTerm(GModel *gm, int iField = 0) :
gmshNodalFemTerm<double>(gm), _iField(iField){}
- virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+ virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
};
#endif
diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
index 5728de1b190602669a89219899e6524541fd0a3a..e6b7870488c73df4d83bc5e3ea31481c690ab6ac 100644
--- a/Numeric/gmshTermOfFormulation.h
+++ b/Numeric/gmshTermOfFormulation.h
@@ -11,7 +11,7 @@
#include <math.h>
#include <map>
#include <vector>
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include "gmshFunction.h"
#include "gmshAssembler.h"
#include "GModel.h"
@@ -48,8 +48,8 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
public:
gmshNodalFemTerm(GModel *gm) : gmshTermOfFormulation<scalar>(gm) {}
virtual ~gmshNodalFemTerm (){}
- virtual void elementMatrix(MElement *e, gmshMatrix<scalar> &m) const = 0;
- virtual void elementVector(MElement *e, gmshVector<scalar> &m) const {
+ virtual void elementMatrix(MElement *e, fullMatrix<scalar> &m) const = 0;
+ virtual void elementVector(MElement *e, fullVector<scalar> &m) const {
m.scale(0.0);
}
void addToMatrix(gmshAssembler<scalar> &lsys) const
@@ -77,7 +77,7 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
{
const int nbR = sizeOfR(e);
const int nbC = sizeOfC(e);
- gmshMatrix<scalar> localMatrix(nbR, nbC);
+ fullMatrix<scalar> localMatrix(nbR, nbC);
elementMatrix(e, localMatrix);
addToMatrix(lsys, localMatrix, e);
}
@@ -86,7 +86,7 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
for (unsigned int i = 0; i < v.size(); i++)
addToMatrix(lsys, v[i]);
}
- void addToMatrix(gmshAssembler<scalar> &lsys, gmshMatrix<scalar> &localMatrix,
+ void addToMatrix(gmshAssembler<scalar> &lsys, fullMatrix<scalar> &localMatrix,
MElement *e) const
{
const int nbR = sizeOfR(e);
@@ -173,7 +173,7 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
MElement *e = ge->getMeshElement (i);
int nbR = sizeOfR(e);
- gmshVector<scalar> V (nbR);
+ fullVector<scalar> V (nbR);
elementVector (e, V);
// assembly
for (int j=0;j<nbR;j++){
diff --git a/Numeric/GmshPredicates.cpp b/Numeric/robustPredicates.cpp
similarity index 99%
rename from Numeric/GmshPredicates.cpp
rename to Numeric/robustPredicates.cpp
index 768f8db47d3d1e4c7521295b0de2e5d7f1d7a576..a861a3094e15308f253965bf9916370534dd96db 100644
--- a/Numeric/GmshPredicates.cpp
+++ b/Numeric/robustPredicates.cpp
@@ -123,7 +123,7 @@
#include <fpu_control.h>
#endif /* LINUX */
-namespace gmsh
+namespace robustPredicates
{
/* On some machines, the exact arithmetic routines might be defeated by the */
@@ -4176,4 +4176,4 @@ REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
return insphereadapt(pa, pb, pc, pd, pe, permanent);
}
-} // namespace gmsh
+} // namespace robustPredicates
diff --git a/Numeric/GmshPredicates.h b/Numeric/robustPredicates.h
similarity index 86%
rename from Numeric/GmshPredicates.h
rename to Numeric/robustPredicates.h
index c4f6535de73b96ec3290e6ea25a04b77fc9f32eb..bee589708f87e92a857fb0fc20b6efb06cc92493 100644
--- a/Numeric/GmshPredicates.h
+++ b/Numeric/robustPredicates.h
@@ -3,11 +3,11 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-#ifndef _GMSH_PREDICATES_H_
-#define _GMSH_PREDICATES_H_
+#ifndef _ROBUST_PREDICATES_H_
+#define _ROBUST_PREDICATES_H_
// namespace necessary to avoid conflicts with predicates used by Tetgen
-namespace gmsh{
+namespace robustPredicates{
double exactinit();
double incircle(double *pa, double *pb, double *pc, double *pd);
double insphere(double *pa, double *pb, double *pc, double *pd, double *pe);
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index d5d22a9da1a54a6fc451ffadf51d98f3d389d01e..4dc7a81523180780af0753cf74418baeba947480 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -340,7 +340,7 @@
#include <time.h>
#include "GmshConfig.h"
#include "GmshMessage.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include "MallocUtils.h"
#include "ListUtils.h"
#include "TreeUtils.h"
@@ -403,7 +403,7 @@ void yyerror(const char *s);
void yymsg(int level, const char *fmt, ...);
void skip_until(const char *skip, const char *until);
int PrintListOfDouble(char *format, List_T *list, char *buffer);
-gmshMatrix<double> ListOfListOfDouble2Matrix(List_T *list);
+fullMatrix<double> ListOfListOfDouble2Matrix(List_T *list);
/* Enabling traces. */
@@ -8728,7 +8728,7 @@ int PrintListOfDouble(char *format, List_T *list, char *buffer)
return 0;
}
-gmshMatrix<double> ListOfListOfDouble2Matrix(List_T *list)
+fullMatrix<double> ListOfListOfDouble2Matrix(List_T *list)
{
int M = List_Nbr(list);
int N = 0;
@@ -8736,7 +8736,7 @@ gmshMatrix<double> ListOfListOfDouble2Matrix(List_T *list)
List_T *line = *(List_T**)List_Pointer_Fast(list, i);
N = std::max(N, List_Nbr(line));
}
- gmshMatrix<double> mat(M, N);
+ fullMatrix<double> mat(M, N);
for(int i = 0; i < M; i++){
List_T *line = *(List_T**)List_Pointer_Fast(list, i);
for(int j = 0; j < List_Nbr(line); j++){
diff --git a/Plugin/CutPlane.cpp b/Plugin/CutPlane.cpp
index 2727bb7dcbec0fd72fb38b1813db1913d98481a2..542f888b7bd29219e18b7561dc3f48ce8b937143 100644
--- a/Plugin/CutPlane.cpp
+++ b/Plugin/CutPlane.cpp
@@ -138,7 +138,7 @@ double GMSH_CutPlanePlugin::levelset(double x, double y, double z, double val) c
CutPlaneOptions_Number[2].def * z + CutPlaneOptions_Number[3].def;
}
-bool GMSH_CutPlanePlugin::geometricalFilter(gmshMatrix<double> *node_positions) const
+bool GMSH_CutPlanePlugin::geometricalFilter(fullMatrix<double> *node_positions) const
{
const double l0 = levelset((*node_positions)(0, 0),
(*node_positions)(0, 1),
diff --git a/Plugin/CutPlane.h b/Plugin/CutPlane.h
index 7aa9efdf6f85f946d3c4827563c88c56e727eec7..c8cf5747dbbd52f4313938584fed4f5b5d1b8578 100644
--- a/Plugin/CutPlane.h
+++ b/Plugin/CutPlane.h
@@ -26,7 +26,7 @@ class GMSH_CutPlanePlugin : public GMSH_LevelsetPlugin
int getNbOptions() const;
StringXNumber *getOption(int iopt);
PView *execute(PView *);
- virtual bool geometricalFilter(gmshMatrix<double> *) const;
+ virtual bool geometricalFilter(fullMatrix<double> *) const;
static double callbackA(int, int, double);
static double callbackB(int, int, double);
diff --git a/Plugin/Plugin.h b/Plugin/Plugin.h
index 43ed85b0a04b0e8635d8f7c346e20f3b58d533a0..d20167ef90740df7988e0a3b3023645bb2c5a88a 100644
--- a/Plugin/Plugin.h
+++ b/Plugin/Plugin.h
@@ -18,7 +18,7 @@
#include "GmshMessage.h"
#include "PView.h"
#include "PViewDataList.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
class PluginDialogBox;
class Vertex;
@@ -88,7 +88,7 @@ class GMSH_PostPlugin : public GMSH_Plugin
// get the data in list format
virtual PViewDataList *getDataList(PView *view, bool showError=true);
virtual void assignSpecificVisibility() const {}
- virtual bool geometricalFilter(gmshMatrix<double> *) const { return true; }
+ virtual bool geometricalFilter(fullMatrix<double> *) const { return true; }
};
// The base class for solver plugins. The idea is to be able to
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index a1cdd26b5f2e6e03431be58131e9c1356fbf6e4a..7367b09b0607c8ab214ff8c67e166a53cf1b665d 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -16,7 +16,7 @@ PViewData::PViewData()
PViewData::~PViewData()
{
if(_adaptive) delete _adaptive;
- for(std::map<int, std::vector<gmshMatrix<double>*> >::iterator it = _interpolation.begin();
+ for(std::map<int, std::vector<fullMatrix<double>*> >::iterator it = _interpolation.begin();
it != _interpolation.end(); it++)
for(unsigned int i = 0; i < it->second.size(); i++)
delete it->second[i];
@@ -73,28 +73,28 @@ void PViewData::setValue(int step, int ent, int ele, int nod, int comp, double v
}
void PViewData::setInterpolationMatrices(int type,
- const gmshMatrix<double> &coefVal,
- const gmshMatrix<double> &expVal)
+ const fullMatrix<double> &coefVal,
+ const fullMatrix<double> &expVal)
{
if(!type || _interpolation[type].size()) return;
- _interpolation[type].push_back(new gmshMatrix<double>(coefVal));
- _interpolation[type].push_back(new gmshMatrix<double>(expVal));
+ _interpolation[type].push_back(new fullMatrix<double>(coefVal));
+ _interpolation[type].push_back(new fullMatrix<double>(expVal));
}
void PViewData::setInterpolationMatrices(int type,
- const gmshMatrix<double> &coefVal,
- const gmshMatrix<double> &expVal,
- const gmshMatrix<double> &coefGeo,
- const gmshMatrix<double> &expGeo)
+ const fullMatrix<double> &coefVal,
+ const fullMatrix<double> &expVal,
+ const fullMatrix<double> &coefGeo,
+ const fullMatrix<double> &expGeo)
{
if(!type || _interpolation[type].size()) return;
- _interpolation[type].push_back(new gmshMatrix<double>(coefVal));
- _interpolation[type].push_back(new gmshMatrix<double>(expVal));
- _interpolation[type].push_back(new gmshMatrix<double>(coefGeo));
- _interpolation[type].push_back(new gmshMatrix<double>(expGeo));
+ _interpolation[type].push_back(new fullMatrix<double>(coefVal));
+ _interpolation[type].push_back(new fullMatrix<double>(expVal));
+ _interpolation[type].push_back(new fullMatrix<double>(coefGeo));
+ _interpolation[type].push_back(new fullMatrix<double>(expGeo));
}
-int PViewData::getInterpolationMatrices(int type, std::vector<gmshMatrix<double>*> &p)
+int PViewData::getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p)
{
if(_interpolation.count(type)){
p = _interpolation[type];
diff --git a/Post/PViewData.h b/Post/PViewData.h
index a80ce3150ff93fc40e2adb76b423eafc5c3b3702..4f04eea04f1df533ae0b3a47b7445bee879c1330 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -10,7 +10,7 @@
#include <vector>
#include <map>
#include "SBoundingBox3d.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#define VAL_INF 1.e200
@@ -34,7 +34,7 @@ class PViewData {
// adaptive visualization data
adaptiveData *_adaptive;
// interpolation matrices, indexed by the type of element
- std::map<int, std::vector<gmshMatrix<double>*> > _interpolation;
+ std::map<int, std::vector<fullMatrix<double>*> > _interpolation;
public:
PViewData();
@@ -180,14 +180,14 @@ class PViewData {
// set/get the interpolation matrices for elements with type "type"
void setInterpolationMatrices(int type,
- const gmshMatrix<double> &coefVal,
- const gmshMatrix<double> &expVal);
+ const fullMatrix<double> &coefVal,
+ const fullMatrix<double> &expVal);
void setInterpolationMatrices(int type,
- const gmshMatrix<double> &coefVal,
- const gmshMatrix<double> &expVal,
- const gmshMatrix<double> &coefGeo,
- const gmshMatrix<double> &expGeo);
- int getInterpolationMatrices(int type, std::vector<gmshMatrix<double>*> &p);
+ const fullMatrix<double> &coefVal,
+ const fullMatrix<double> &expVal,
+ const fullMatrix<double> &coefGeo,
+ const fullMatrix<double> &expGeo);
+ int getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p);
inline bool haveInterpolationMatrices(){ return !_interpolation.empty(); }
// smooth the data in the view (makes it C0)
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index bf951be60ce84fcd5bbb043c69e2c183bfa9bcbe..1a9f11a5f647a65f7390ef4ee05975541df13734 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -6,7 +6,7 @@
#include "PViewDataList.h"
#include "GmshMessage.h"
#include "GmshDefines.h"
-#include "FunctionSpace.h"
+#include "functionSpace.h"
#include "Numeric.h"
#include "SmoothData.h"
#include "Context.h"
@@ -148,7 +148,7 @@ void PViewDataList::_stat(std::vector<double> &list, int nbcomp, int nbelm,
int nbval = nbcomp * nbnod;
if(haveInterpolationMatrices()){
- std::vector<gmshMatrix<double>*> im;
+ std::vector<fullMatrix<double>*> im;
int nim = getInterpolationMatrices(type, im);
if(nim == 4)
nbnod = im[2]->size1();
@@ -201,7 +201,7 @@ void PViewDataList::_setLast(int ele, int dim, int nbnod, int nbcomp, int nbedg,
std::vector<double> &list, int nblist)
{
if(haveInterpolationMatrices()){
- std::vector<gmshMatrix<double>*> im;
+ std::vector<fullMatrix<double>*> im;
if(getInterpolationMatrices(type, im) == 4)
nbnod = im[2]->size1();
}
@@ -529,11 +529,11 @@ bool PViewDataList::combineSpace(nameData &nd)
}
// copy interpolation marices
- for(std::map<int, std::vector<gmshMatrix<double>*> >::iterator it =
+ for(std::map<int, std::vector<fullMatrix<double>*> >::iterator it =
l->_interpolation.begin(); it != l->_interpolation.end(); it++)
if(_interpolation[it->first].empty())
for(unsigned int i = 0; i < it->second.size(); i++)
- _interpolation[it->first].push_back(new gmshMatrix<double>(*it->second[i]));
+ _interpolation[it->first].push_back(new fullMatrix<double>(*it->second[i]));
// merge elememts
dVecMerge(l->SP, SP); NbSP += l->NbSP; dVecMerge(l->VP, VP); NbVP += l->NbVP;
@@ -624,11 +624,11 @@ bool PViewDataList::combineTime(nameData &nd)
}
NbT2 = data[0]->NbT2;
NbT3 = data[0]->NbT3;
- for(std::map<int, std::vector<gmshMatrix<double>*> >::iterator it =
+ for(std::map<int, std::vector<fullMatrix<double>*> >::iterator it =
data[0]->_interpolation.begin(); it != data[0]->_interpolation.end(); it++)
if(_interpolation[it->first].empty())
for(unsigned int i = 0; i < it->second.size(); i++)
- _interpolation[it->first].push_back(new gmshMatrix<double>(*it->second[i]));
+ _interpolation[it->first].push_back(new fullMatrix<double>(*it->second[i]));
// merge values for all element types
for(int i = 0; i < 24; i++){
@@ -784,7 +784,7 @@ void PViewDataList::setOrder2(int type)
case TYPE_PRI: typeMSH = MSH_PRI_18; break;
case TYPE_PYR: typeMSH = MSH_PYR_14; break;
}
- const gmshFunctionSpace *fs = &gmshFunctionSpaces::find(typeMSH);
+ const functionSpace *fs = &functionSpaces::find(typeMSH);
if(!fs){
Msg::Error("Could not find function space for element type %d", typeMSH);
return;
diff --git a/Post/adaptiveData.cpp b/Post/adaptiveData.cpp
index a5d7b9c9c05c7948f47777f8b6741759f7f90881..0571f954210d00ec30bd9f671a1ece4044d15d73 100644
--- a/Post/adaptiveData.cpp
+++ b/Post/adaptiveData.cpp
@@ -50,9 +50,9 @@ static void cleanElement()
T::allPoints.clear();
}
-static void computeShapeFunctions(gmshMatrix<double> *coeffs, gmshMatrix<double> *eexps,
- double u, double v, double w, gmshVector<double> *sf,
- gmshVector<double> *tmp)
+static void computeShapeFunctions(fullMatrix<double> *coeffs, fullMatrix<double> *eexps,
+ double u, double v, double w, fullVector<double> *sf,
+ fullVector<double> *tmp)
{
for(int i = 0; i < eexps->size1(); i++) {
(*tmp)(i) = pow(u, (*eexps)(i, 0));
@@ -830,7 +830,7 @@ void adaptivePrism::recurError(adaptivePrism *p, double AVG, double tol)
}
template <class T>
-adaptiveElements<T>::adaptiveElements(std::vector<gmshMatrix<double>*> &p)
+adaptiveElements<T>::adaptiveElements(std::vector<fullMatrix<double>*> &p)
: _coeffsVal(0), _eexpsVal(0), _interpolVal(0),
_coeffsGeom(0), _eexpsGeom(0), _interpolGeom(0)
{
@@ -864,15 +864,15 @@ void adaptiveElements<T>::init(int level)
int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
if(_interpolVal) delete _interpolVal;
- _interpolVal = new gmshMatrix<double>(T::allPoints.size(), numVals);
+ _interpolVal = new fullMatrix<double>(T::allPoints.size(), numVals);
if(_interpolGeom) delete _interpolGeom;
- _interpolGeom = new gmshMatrix<double>(T::allPoints.size(), numNodes);
+ _interpolGeom = new fullMatrix<double>(T::allPoints.size(), numNodes);
- gmshVector<double> sfv(numVals), *tmpv = 0;
- gmshVector<double> sfg(numNodes), *tmpg = 0;
- if(_eexpsVal) tmpv = new gmshVector<double>(_eexpsVal->size1());
- if(_eexpsGeom) tmpg = new gmshVector<double>(_eexpsGeom->size1());
+ fullVector<double> sfv(numVals), *tmpv = 0;
+ fullVector<double> sfg(numNodes), *tmpg = 0;
+ if(_eexpsVal) tmpv = new fullVector<double>(_eexpsVal->size1());
+ if(_eexpsGeom) tmpg = new fullVector<double>(_eexpsGeom->size1());
int i = 0;
for(std::set<adaptivePoint>::iterator it = T::allPoints.begin();
@@ -937,7 +937,7 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
double t1 = GetTimeInSeconds();
#endif
- gmshVector<double> val(numVals), res(numPoints);
+ fullVector<double> val(numVals), res(numPoints);
if(numComp == 1){
for(int i = 0; i < numVals; i++)
val(i) = values[i].v[0];
@@ -957,10 +957,10 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
}
if(onlyComputeMinMax) return;
- gmshMatrix<double> *resxyz = 0;
+ fullMatrix<double> *resxyz = 0;
if(numComp == 3){
- gmshMatrix<double> valxyz(numVals, 3);
- resxyz = new gmshMatrix<double>(numPoints, 3);
+ fullMatrix<double> valxyz(numVals, 3);
+ resxyz = new fullMatrix<double>(numPoints, 3);
for(int i = 0; i < numVals; i++){
valxyz(i, 0) = values[i].v[0];
valxyz(i, 1) = values[i].v[1];
@@ -976,7 +976,7 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
return;
}
- gmshMatrix<double> xyz(numNodes, 3), XYZ(numPoints, 3);
+ fullMatrix<double> xyz(numNodes, 3), XYZ(numPoints, 3);
for(int i = 0; i < numNodes; i++){
xyz(i, 0) = coords[i].c[0];
xyz(i, 1) = coords[i].c[1];
@@ -1134,7 +1134,7 @@ adaptiveData::adaptiveData(PViewData *data)
_tetrahedra(0), _hexahedra(0), _prisms(0)
{
_outData = new PViewDataList();
- std::vector<gmshMatrix<double>*> p;
+ std::vector<fullMatrix<double>*> p;
if(_inData->getNumLines()){
_inData->getInterpolationMatrices(TYPE_LIN, p);
_lines = new adaptiveElements<adaptiveLine>(p);
diff --git a/Post/adaptiveData.h b/Post/adaptiveData.h
index 4f7550624559807cb9aae7adabf418651ace4073..aba96a5a45a3fe16b12c3973503679f912b05e9a 100644
--- a/Post/adaptiveData.h
+++ b/Post/adaptiveData.h
@@ -9,7 +9,7 @@
#include <list>
#include <set>
#include <vector>
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
class PViewData;
class PViewDataList;
@@ -53,7 +53,7 @@ class adaptiveLine {
{
return (p[0]->val + p[1]->val) / 2.;
}
- inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
+ inline static void GSF(double u, double v, double w, fullVector<double> &sf)
{
sf(0) = (1 - u) / 2.;
sf(1) = (1 + u) / 2.;
@@ -85,7 +85,7 @@ class adaptiveTriangle {
{
return (p[0]->val + p[1]->val + p[2]->val) / 3.;
}
- inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
+ inline static void GSF(double u, double v, double w, fullVector<double> &sf)
{
sf(0) = 1. - u - v;
sf(1) = u;
@@ -120,7 +120,7 @@ class adaptiveQuadrangle {
{
return (p[0]->val + p[1]->val + p[2]->val + p[3]->val) / 4.;
}
- inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
+ inline static void GSF(double u, double v, double w, fullVector<double> &sf)
{
sf(0) = 0.25 * (1. - u) * (1. - v);
sf(1) = 0.25 * (1. + u) * (1. - v);
@@ -161,7 +161,7 @@ class adaptivePrism {
{
return (p[0]->val + p[1]->val + p[2]->val + p[3]->val + p[4]->val + p[5]->val) / 6.;
}
- inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
+ inline static void GSF(double u, double v, double w, fullVector<double> &sf)
{
sf(0) = (1. - u - v) * (1 - w) / 2;
sf(1) = u * (1-w)/2;
@@ -200,7 +200,7 @@ class adaptiveTetrahedron {
{
return (p[0]->val + p[1]->val + p[2]->val + p[3]->val) / 4.;
}
- inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
+ inline static void GSF(double u, double v, double w, fullVector<double> &sf)
{
sf(0) = 1. - u - v - w;
sf(1) = u;
@@ -243,7 +243,7 @@ class adaptiveHexahedron {
return (p[0]->val + p[1]->val + p[2]->val+ p[3]->val +
p[4]->val + p[5]->val + p[6]->val+ p[7]->val) / 8.;
}
- inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
+ inline static void GSF(double u, double v, double w, fullVector<double> &sf)
{
sf(0) = 0.125 * (1 - u) * (1 - v) * (1 - w);
sf(1) = 0.125 * (1 + u) * (1 - v) * (1 - w);
@@ -285,10 +285,10 @@ class PValues{
template <class T>
class adaptiveElements {
private:
- gmshMatrix<double> *_coeffsVal, *_eexpsVal, *_interpolVal;
- gmshMatrix<double> *_coeffsGeom, *_eexpsGeom, *_interpolGeom;
+ fullMatrix<double> *_coeffsVal, *_eexpsVal, *_interpolVal;
+ fullMatrix<double> *_coeffsGeom, *_eexpsGeom, *_interpolGeom;
public:
- adaptiveElements(std::vector<gmshMatrix<double>*> &interpolationMatrices);
+ adaptiveElements(std::vector<fullMatrix<double>*> &interpolationMatrices);
~adaptiveElements();
// create the _interpolVal and _interpolGeom matrices at the given
// refinement level
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index 06eec04a06f4e916faa1d0ed2664e4aeede9e0dd..44e7293cee7f33f70692ccd8076cdce7b03fa903 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -6,7 +6,7 @@
#include "elasticityTerm.h"
#include "Numeric.h"
-void elasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+void elasticityTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
@@ -31,10 +31,10 @@ void elasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
{ 0, 0, 0, 0, C44, 0},
{ 0, 0, 0, 0, 0, C44} };
- gmshMatrix<double> H(6, 6);
- gmshMatrix<double> B(6, 3 * nbNodes);
- gmshMatrix<double> BTH(3 * nbNodes, 6);
- gmshMatrix<double> BT(3 * nbNodes, 6);
+ fullMatrix<double> H(6, 6);
+ fullMatrix<double> B(6, 3 * nbNodes);
+ fullMatrix<double> BTH(3 * nbNodes, 6);
+ fullMatrix<double> BT(3 * nbNodes, 6);
for (int i = 0; i < 6; i++)
for (int j = 0; j < 6; j++)
H(i, j) = C[i][j];
@@ -75,7 +75,7 @@ void elasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
}
}
-void elasticityTerm::elementVector(MElement *e, gmshVector<double> &m) const
+void elasticityTerm::elementVector(MElement *e, fullVector<double> &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
diff --git a/Solver/elasticityTerm.h b/Solver/elasticityTerm.h
index db1c4574586dbe11d421dc7b6912778af36639c2..e054c04fa3a2b542b0a0aed0269b53d8928340cb 100644
--- a/Solver/elasticityTerm.h
+++ b/Solver/elasticityTerm.h
@@ -10,7 +10,7 @@
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
class elasticityTerm : public femTerm<double,double> {
protected:
@@ -34,8 +34,8 @@ class elasticityTerm : public femTerm<double,double> {
elasticityTerm(GModel *gm, double E, double nu, int iField) :
femTerm<double,double>(gm), _E(E), _nu(nu), _iField(iField){}
void setVector(const SVector3 &f) {_volumeForce = f;}
- void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
- void elementVector(MElement *e, gmshVector<double> &m) const;
+ void elementMatrix(MElement *e, fullMatrix<double> &m) const;
+ void elementVector(MElement *e, fullVector<double> &m) const;
};
#endif
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index 001109e203ef1fa96f1382ee4fd3a6bf303db7f7..51327f07048f1b0718764c079edc709fe994344c 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -9,7 +9,7 @@
#include <math.h>
#include <map>
#include <vector>
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include "gmshFunction.h"
#include "dofManager.h"
#include "GModel.h"
@@ -37,8 +37,8 @@ class femTerm {
public:
femTerm(GModel *gm) : _gm(gm) {}
virtual ~femTerm (){}
- virtual void elementMatrix(MElement *e, gmshMatrix<dataMat> &m) const = 0;
- virtual void elementVector(MElement *e, gmshVector<dataVec> &m) const
+ virtual void elementMatrix(MElement *e, fullMatrix<dataMat> &m) const = 0;
+ virtual void elementVector(MElement *e, fullVector<dataVec> &m) const
{
m.scale(0.0);
}
@@ -53,12 +53,12 @@ class femTerm {
{
const int nbR = sizeOfR(e);
const int nbC = sizeOfC(e);
- gmshMatrix<dataMat> localMatrix(nbR, nbC);
+ fullMatrix<dataMat> localMatrix(nbR, nbC);
elementMatrix(e, localMatrix);
addToMatrix(dm, localMatrix, e);
}
void addToMatrix(dofManager<dataVec,dataMat> &dm,
- gmshMatrix<dataMat> &localMatrix,
+ fullMatrix<dataMat> &localMatrix,
MElement *e) const
{
const int nbR = localMatrix.size1();
@@ -122,7 +122,7 @@ class femTerm {
for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
MElement *e = ge->getMeshElement (i);
int nbR = sizeOfR(e);
- gmshVector<dataVec> V (nbR);
+ fullVector<dataVec> V (nbR);
elementVector (e, V);
// assembly
for (int j=0;j<nbR;j++)dm.assemble(getLocalDofR(e,j),V(j));
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index e88cfe3e8e8ad4f6dc2f7fc52747bedc23e255af..73a4dc159a2c843ddfd3531e9c7c2dcd6d8dcaca 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -12,7 +12,7 @@
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include "Numeric.h"
// \nabla \cdot k \nabla U + a U
@@ -39,7 +39,7 @@ class helmoltzTerm : public femTerm<scalar,scalar> {
gmshFunction<scalar> *k,
gmshFunction<scalar> *a) :
femTerm<scalar,scalar>(gm), _k(k), _a(a), _iFieldR(iFieldR), _iFieldC(iFieldC){}
- virtual void elementMatrix(MElement *e, gmshMatrix<scalar> &m) const
+ virtual void elementMatrix(MElement *e, fullMatrix<scalar> &m) const
{
// compute integration rule
const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() :
diff --git a/Solver/linearSystemFull.h b/Solver/linearSystemFull.h
index 831fc60110485d8ea3910c9f3ba0ebdc811594de..932af4443941c0c086ea00b1e8149254b5b09ace 100644
--- a/Solver/linearSystemFull.h
+++ b/Solver/linearSystemFull.h
@@ -10,24 +10,24 @@
#include "GmshMessage.h"
#include "linearSystem.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include <stdlib.h>
#include <set>
template <class scalar>
class linearSystemFull : public linearSystem<scalar> {
private:
- gmshMatrix<scalar> *_a;
- gmshVector<scalar> *_b, *_x;
+ fullMatrix<scalar> *_a;
+ fullVector<scalar> *_b, *_x;
public :
linearSystemFull() : _a(0), _b(0), _x(0){}
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int _nbRows)
{
clear();
- _a = new gmshMatrix<scalar>(_nbRows, _nbRows);
- _b = new gmshVector<scalar>(_nbRows);
- _x = new gmshVector<scalar>(_nbRows);
+ _a = new fullMatrix<scalar>(_nbRows, _nbRows);
+ _b = new fullVector<scalar>(_nbRows);
+ _x = new fullVector<scalar>(_nbRows);
}
virtual ~linearSystemFull()
{
diff --git a/utils/api_demos/mainLaplace.cpp b/utils/api_demos/mainLaplace.cpp
deleted file mode 100644
index c269c129d5ed62fa61291a12bebbc711fe974255..0000000000000000000000000000000000000000
--- a/utils/api_demos/mainLaplace.cpp
+++ /dev/null
@@ -1,877 +0,0 @@
-#include <stdio.h>
-#include <gmsh/Gmsh.h>
-#include <gmsh/GModel.h>
-#include <gmsh/MElement.h>
-#include <gmsh/GmshMatrix.h>
-class gLevelset;
-#ifdef HAVE_GLEVELSETS
-#include <Levelset.h>
-#include <Integration3D.h>
-#endif
-
-// My First test with Gmsh's API
-// A Laplace Equation (booh)
-// Vector ort scalar, in any dimension and with
-// any polynomial order
-
-// This code will quantify the gain of assembling matrices
-// either in reference coordinate system, or in the real one
-
-typedef Double_Matrix gmshSmallMatrix ;
-#ifdef HAVE_SPARSKIT
-#include "CSR_Matrix.h"
-#include "CSR_Vector.h"
-typedef CSR_Matrix gmshMatrix ;
-typedef CSR_Vector gmshVector ;
-inline double getMatrix ( int i, int j, gmshMatrix &m){
- return m.GetMatrix(i+1,j+1);
-}
-inline void addMatrix ( int i, int j, double val, gmshMatrix &m){
- if (val != 0.0)m.AddMatrix(i+1,j+1,val);
-}
-inline void addVector ( int i, double val, gmshVector &v){
- if (val != 0.0)v.AddVal(i+1,val);
-}
-inline void zeroVector ( gmshVector &v){
- v.ZeroArray();
-}
-inline double getVector ( int i, gmshVector &v){
- return v.GetVal(i+1);
-}
-inline void SystemSolve ( gmshMatrix &a, gmshVector &x, gmshVector &b){
- a.EndOfAssembly();
- SPARSKIT_LINEAR_SOLVER_ ( "rcmk","ilut","gmres",a,b,x);
-}
-#else
-typedef Double_Matrix gmshMatrix ;
-typedef Double_Vector gmshVector ;
-inline double getMatrix ( int i, int j, gmshMatrix &m){
- return m(i,j);
-}
-inline void addMatrix ( int i, int j, double val, gmshMatrix &m){
- m(i,j) += val;
-}
-inline void addVector ( int i, double val, gmshVector &v){
- v(i) += val;
-}
-inline double getVector ( int i, gmshVector &v){
- return v(i);
-}
-inline void SystemSolve ( gmshMatrix &a, gmshVector &x, gmshVector &b){
- a.lu_solve (b,x);
-}
-inline void zeroVector ( gmshVector &v){
- v.set_all(0.0);
-}
-
-#endif
-
-extern double inv3x3 (double a[3][3], double a[3][3]);
-
-// An ad hoc structure to store the triangulation
-// issued of a levelset cut
-// num is the new vertex
-// num1 and num2 are vertices of the cut edge in the initial mesh
-class gmshTopoVertex
-{
- static double tol;
- int num;
- double x,y,z;
- MVertex *v;
-public:
- gmshTopoVertex (double _x, double _y, double _z):
- num(-1),x(_x),y(_y),z(_z)
- {
- }
- void setMVertex(MVertex *_v){v=_v;}
- MVertex* getMVertex()const {return v;}
- void setNum(int n){num=n;}
- int getNum()const {return num;}
- bool operator < (const gmshTopoVertex & other) const{
- if (x - other.x < -tol) return true;
- if (x - other.x > tol) return false;
- if (y - other.y < -tol) return true;
- if (y - other.y > tol) return false;
- if (z - other.z < -tol) return true;
- return false;
- }
- double X() const {return x;};
- double Y() const {return y;};
- double Z() const {return z;};
-};
-
-double gmshTopoVertex::tol = 1.e-9;
-
-// try to add a "poisson" term INSIDE the elements
-#ifdef HAVE_GLEVELSETS
-void IntegrationPoints (MElement*e,
- gLevelset*ls,
- int integrationOrder,
- std::vector<IntegrationPoint>&ipV ,
- std::vector<IntegrationPoint>&ipS )
-{
- switch ( e->getTypeForMSH()){
- case MSH_TET_4 :
- {
- vector<CuttingPoint> cp;
- vector<Tetra> subTetras;
- vector<Quad> surfQuads;
- vector<Triangle> surfTriangles;
- Tetra T (e->getVertex(0)->x(),
- e->getVertex(0)->y(),
- e->getVertex(0)->z(),
- e->getVertex(1)->x(),
- e->getVertex(1)->y(),
- e->getVertex(1)->z(),
- e->getVertex(2)->x(),
- e->getVertex(2)->y(),
- e->getVertex(2)->z(),
- e->getVertex(3)->x(),
- e->getVertex(3)->y(),
- e->getVertex(3)->z());
- T.cut (*ls, ipV, ipS, cp, integrationOrder,integrationOrder,integrationOrder, subTetras, surfQuads, surfTriangles);
- }
- break;
-
- default :
- throw;
- }
-}
-void gmshCutMesh (MElement*e,
- gLevelset*ls,
- std::set<gmshTopoVertex> &tVerts,
- std::multimap<MElement*,MElement*> &cut)
-{
- switch ( e->getTypeForMSH()){
- case MSH_TET_4 :
- {
- vector<CuttingPoint> cp;
- vector<Tetra> subTetras;
- vector<Quad> surfQuads;
- vector<Triangle> surfTriangles;
- int integrationOrder = 1;
- std::vector<IntegrationPoint> ipV;
- std::vector<IntegrationPoint> ipS;
-
- Tetra T (e->getVertex(0)->x(),
- e->getVertex(0)->y(),
- e->getVertex(0)->z(),
- e->getVertex(1)->x(),
- e->getVertex(1)->y(),
- e->getVertex(1)->z(),
- e->getVertex(2)->x(),
- e->getVertex(2)->y(),
- e->getVertex(2)->z(),
- e->getVertex(3)->x(),
- e->getVertex(3)->y(),
- e->getVertex(3)->z());
- T.cut (*ls, ipV, ipS, cp, integrationOrder,integrationOrder,integrationOrder, subTetras, surfQuads, surfTriangles);
- for (int i=0;i<surfTriangles.size();i++){
- gmshTopoVertex v0 ( surfTriangles[i].x(0),surfTriangles[i].y(0),surfTriangles[i].z(0));
- gmshTopoVertex v1 ( surfTriangles[i].x(1),surfTriangles[i].y(1),surfTriangles[i].z(1));
- gmshTopoVertex v2 ( surfTriangles[i].x(2),surfTriangles[i].y(2),surfTriangles[i].z(2));
- std::set<gmshTopoVertex>::const_iterator it0 = tVerts.find(v0);
- std::set<gmshTopoVertex>::iterator it1 = tVerts.find(v1);
- std::set<gmshTopoVertex>::iterator it2 = tVerts.find(v2);
- MVertex *mv0,*mv1,*mv2;
- if (it0 == tVerts.end()){
- v0.setNum(tVerts.size());
- mv0 = new MVertex (v0.X(),v0.Y(),v0.Z(),0,v0.getNum());
- v0.setMVertex (mv0);
- tVerts.insert(v0);
- }
- else mv0 = it0->getMVertex();
- if (it1 == tVerts.end()){
- v1.setNum(tVerts.size());
- mv1 = new MVertex (v1.X(),v1.Y(),v1.Z(),0,v1.getNum());
- v1.setMVertex (mv1);
- tVerts.insert(v1);
- }
- else mv1 = it1->getMVertex();
- if (it2 == tVerts.end()){
- v2.setNum(tVerts.size());
- mv2 = new MVertex (v2.X(),v2.Y(),v2.Z(),0,v2.getNum());
- v2.setMVertex (mv2);
- tVerts.insert(v2);
- }
- else mv2 = it2->getMVertex();
-
- MTriangle *tri = new MTriangle(mv0,mv1,mv2);
- cut.insert(std::make_pair(e,tri));
- }
- }
- break;
-
- default :
- throw;
- }
-}
-
-
-void gmshBuildSurfaceMesh(GModel *_gm,
- gLevelset *ls,
- std::set<gmshTopoVertex> &tVerts,
- std::multimap<MElement*,MElement*> &cut) {
- if (_gm->getNumRegions()){
- for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
- for (int i = 0;i<(*it)->getNumMeshElements();i++){
- MElement *e = (*it)->getMeshElement(i);
- gmshCutMesh (e,ls,tVerts,cut);
- }
- }
- }
- else throw;
-}
-
-#endif
-
-
-class gmshVertexEvaluator
-{
- double _val;
- public :
- gmshVertexEvaluator(double val = 0):_val(val){}
- virtual double operator () (MVertex *) const {return _val;}
- virtual double operator () (double x, double y, double z) const {return _val;}
-};
-
-class gmshTermOfFormulation {
-protected:
- GModel *_gm;
-public:
- gmshTermOfFormulation (GModel *gm) : _gm(gm){}
- virtual void addToRightHandSide (gmshVector &r) const = 0;
- virtual void addToJacobian (gmshMatrix &J) const = 0;
-};
-
-class gmshNodalFemTerm : public gmshTermOfFormulation {
-protected:
- virtual int sizeOfC(MElement*) const = 0;
- virtual int sizeOfR(MElement*) const = 0;
- virtual int localToGlobalR (MElement *e, int iVertex, int icomp) const = 0;
- virtual int localToGlobalR (MElement *e, int i) const = 0;
- virtual int localToGlobalC (MElement *e, int j) const = 0;
- std::map<std::pair<int,int>, double> dirichlet;
- gLevelset *_ls;
-public:
- virtual void elementMatrix ( MElement *e, gmshSmallMatrix & m) const = 0;
- virtual void elementVector ( MElement *e, gmshVector & m) const = 0;
- void addToJacobian (gmshMatrix &J, GEntity *ge) const;
- void addToRightHandSide (gmshVector &J, GEntity *ge) const;
- gmshNodalFemTerm (GModel *gm) : gmshTermOfFormulation(gm),_ls(0){}
- void addToRightHandSide (gmshVector &r) const;
- void addToJacobian (gmshMatrix &J) const;
- void addToJacobian (gmshMatrix &Jac, gmshSmallMatrix &localMatrix, MElement *e) const;
- void uglyDirichlet (gmshMatrix &J,gmshVector &r) const;
- void addDirichlet (int physical, int dim, int comp, const gmshVertexEvaluator & e);
- void addNeumann (int physical, int dim, int icomp, const gmshVertexEvaluator & e, gmshVector &r);
- void setLevelset(gLevelset *ls) {_ls = ls;}
- virtual int sizeOfR () const = 0;
- virtual int sizeOfC () const = 0;
- virtual void getComponents (int iVertex, int *components) const = 0;
-};
-
-void gmshNodalFemTerm::addDirichlet(int physical, int dim, int comp, const gmshVertexEvaluator & e) {
- std::vector<MVertex *> v;
- _gm->getMeshVertices (physical, dim,v);
- for (int i=0;i<v.size();i++)
- dirichlet[std::make_pair(v[i]->getNum(),comp)] = e(v[i]);
-}
-
-void gmshNodalFemTerm::addNeumann(int physical, int dim, int comp, const gmshVertexEvaluator & fct, gmshVector &r) {
- std::map<int, std::vector<GEntity*> > groups[4];
- _gm->getPhysicalGroups(groups);
- std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical);
- if (it == groups[dim].end())return;
- double jac[3][3];
-
- for (int i=0;i<it->second.size();++i){
- GEntity *ge = it->second[i];
- for (int j=0; j<ge->getNumMeshElements();j++){
- MElement *e = ge->getMeshElement(j);
- int integrationOrder = 2*e->getPolynomialOrder();
- int nbNodes = e->getNumVertices();
- int npts;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- for (int ip=0;ip<npts;ip++){
- const double u = GP[ip].pt[0];
- const double v = GP[ip].pt[1];
- const double w = GP[ip].pt[2];
- const double weight = GP[ip].weight;
- const double detJ = e->getJacobian (u,v,w,jac);
- double x = 0;
- double y = 0;
- double z = 0;
- for (int k=0;k<nbNodes;k++){
- double sf;
- e->getShapeFunction (k,u,v,w,sf);
- x += e->getVertex(k)->x() * sf;
- y += e->getVertex(k)->y() * sf;
- z += e->getVertex(k)->z() * sf;
- }
- const double FCT = fct (x,y,z);
- for (int k=0;k<nbNodes;k++){
- double sf;
- e->getShapeFunction (k,u,v,w,sf);
- // printf("l2g (%d,%d) = %d (%d)\n", k,comp,localToGlobalR(e,k,comp),r.size());
- addVector(localToGlobalR(e,k,comp),detJ * weight * sf * FCT,r);
- }
- }
- }
- }
-}
-
-void gmshNodalFemTerm::uglyDirichlet (gmshMatrix &J,gmshVector &r) const{
- std::map<std::pair<int,int>, double>::const_iterator it = dirichlet.begin();
- int C[256];
- for ( ; it != dirichlet.end() ; ++it){
- getComponents (it->first.first,C);
- double val = it->second;
- int row = C[it->first.second];
- // printf("vertex %d comp %d iRow = %d\n",it->first.first,it->first.second,row);
- const double BIG = getMatrix(row,row,J) * 1.e8;
- addMatrix(row,row,BIG,J);
- addVector(row,BIG * val,r);
- }
-}
-
-// first symmetric implementation, will be better next
-void gmshNodalFemTerm::addToJacobian (gmshMatrix &J) const {
- if (_gm->getNumRegions()){
- for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
- addToJacobian(J,*it);
- }
- }
- else if(_gm->getNumFaces()){
- for(GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it){
- addToJacobian(J,*it);
- }
- }
-}
-
-// first symmetric implementation, will be better next
-void gmshNodalFemTerm::addToRightHandSide (gmshVector &RHS) const {
- if (_gm->getNumRegions()){
- for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
- addToRightHandSide(RHS,*it);
- }
- }
- else if(_gm->getNumFaces()){
- for(GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it){
- addToRightHandSide(RHS,*it);
- }
- }
-}
-
-void gmshNodalFemTerm::addToRightHandSide (gmshVector &RHS, GEntity *ge) const {
- for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement (i);
- int nbR = sizeOfR(e);
- gmshVector V (nbR);
- elementVector (e, V);
- // assembly
- for (int j=0;j<nbR;j++){
- addVector(localToGlobalR (e,j),getVector(j,V),RHS);
- }
- }
-}
-
-void gmshNodalFemTerm::addToJacobian (gmshMatrix &Jac, gmshSmallMatrix &localMatrix, MElement *e) const {
- const int nbR = sizeOfR(e);
- const int nbC = sizeOfC(e);
- for (int j=0;j<nbR;j++){
- int J = localToGlobalR (e,j);
- for (int k=0;k<nbC;k++){
- int K = localToGlobalC (e,k);
- addMatrix(J,K,localMatrix (j,k),Jac);
- }
- }
-}
-
-void gmshNodalFemTerm::addToJacobian (gmshMatrix &Jac, GEntity *ge) const {
- for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement (i);
- const int nbR = sizeOfR(e);
- const int nbC = sizeOfC(e);
- gmshSmallMatrix localMatrix (nbR,nbC);
- elementMatrix (e,localMatrix);
- // assembly
- for (int j=0;j<nbR;j++){
- int J = localToGlobalR (e,j);
- for (int k=0;k<nbC;k++){
- int K = localToGlobalC (e,k);
- addMatrix(J,K,localMatrix (j,k),Jac);
- }
- }
- }
-}
-// LAPLACE
-class gmshLaplaceTerm : public gmshNodalFemTerm {
- const double _diffusivity;
-protected:
- virtual int sizeOfR (MElement *e) const {return e->getNumVertices();}
- virtual int sizeOfC (MElement *e) const {return e->getNumVertices();}
- virtual int localToGlobalR (MElement *e, int i) const {return e->getVertex(i)->getNum()-1;}
- virtual int localToGlobalC (MElement *e, int j) const {return e->getVertex(j)->getNum()-1;}
- virtual int localToGlobalR (MElement *e, int iVertex, int icomp) const {return localToGlobalR (e, iVertex);}
-public:
- gmshLaplaceTerm (GModel *gm, double diffusivity = 1.0) : gmshNodalFemTerm(gm),_diffusivity (diffusivity){}
- void elementMatrix ( MElement *e, gmshSmallMatrix & m) const;
- void elementVector ( MElement *e, gmshVector & v) const;
- virtual int sizeOfR () const {return _gm->getNumMeshVertices();}
- virtual int sizeOfC () const {return _gm->getNumMeshVertices();}
- virtual void getComponents (int iVertex, int *components) const {components[0] = iVertex - 1;};
-};
-
-
-void gmshLaplaceTerm:: elementMatrix ( MElement *e, gmshSmallMatrix & m) const{
- int nbNodes = e->getNumVertices();
- int integrationOrder = 2*(e->getPolynomialOrder()-1);
- int npts;
- IntPt *GP;
- double jac [3][3];
- double invjac [3][3];
- double Grads[256][3],grad[3];
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
-
- m.set_all(0.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);
- inv3x3 ( jac, invjac) ;
- for (int j=0;j<nbNodes;j++){
- e->getGradShapeFunction (j,u,v,w,grad);
- Grads[j][0] = invjac[0][0] * grad[0] + invjac[0][1] * grad[1] + invjac[0][2] * grad[2];
- Grads[j][1] = invjac[1][0] * grad[0] + invjac[1][1] * grad[1] + invjac[1][2] * grad[2];
- Grads[j][2] = invjac[2][0] * grad[0] + invjac[2][1] * grad[1] + invjac[2][2] * grad[2];
- }
- for (int j=0;j<nbNodes;j++){
- for (int k=0;k<=j;k++){
- m(j,k) += (Grads[j][0]*Grads[k][0]+
- Grads[j][1]*Grads[k][1]+
- Grads[j][2]*Grads[k][2]) * weight * detJ * _diffusivity;
- }
- }
- }
- for (int j=0;j<nbNodes;j++)
- for (int k=0;k<j;k++)
- m(k,j) = m(j,k);
-}
-
-
-
-void gmshLaplaceTerm:: elementVector ( MElement *e, gmshVector & m) const{
- int nbNodes = e->getNumVertices();
- int integrationOrder = 2*(e->getPolynomialOrder()-1);
- double jac [3][3];
- zeroVector(m);
-#ifdef HAVE_GLEVELSETS
- std::vector<IntegrationPoint> ipV,ipS;
- IntegrationPoints (e,_ls,integrationOrder,ipV,ipS);
- for (int i=0;i<ipS.size();i++){
- const double u = ipS[i].xl();
- const double v = ipS[i].yl();
- const double w = ipS[i].zl();
- const double weight = ipS[i].weight();
- // const double detJ = e->getJacobian (u,v,w,jac);
- for (int j=0;j<nbNodes;j++){
- double sf ; e->getShapeFunction (j,u,v,w,sf);
- addVector(j,weight*sf,m);
- }
- }
-#else
- int npts;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- 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);
- for (int j=0;j<nbNodes;j++){
- double sf ; e->getShapeFunction (j,u,v,w,sf);
- addVector(j,weight*detJ*sf,m);
- }
- }
-#endif
-}
-
-// LAPLACE
-class gmshLaplaceVectorTerm : public gmshNodalFemTerm {
- double _E,_nu;
-protected:
- virtual int sizeOfR (MElement *e) const {return 3*e->getNumVertices();}
- virtual int sizeOfC (MElement *e) const {return 3*e->getNumVertices();}
- virtual int localToGlobalR (MElement *e, int i) const {
- int iComp = i/e->getNumVertices();
- int ithLocalVertex = i%e->getNumVertices();
- // printf("%d -> %d %d\n",i,ithLocalVertex,iComp);
- return 3*(e->getVertex(ithLocalVertex)->getNum()-1)+iComp;
- }
- virtual int localToGlobalR (MElement *e, int iVertex, int icomp) const {
- return localToGlobalR (e, icomp*e->getNumVertices()+iVertex);
- }
- virtual int localToGlobalC (MElement *e, int j) const {return localToGlobalR (e,j);}
-public:
- gmshLaplaceVectorTerm (GModel *gm, double E, double nu) : gmshNodalFemTerm(gm),_E(E),_nu(nu){}
- void elementMatrix ( MElement *e, gmshSmallMatrix & m) const;
- void elementVector ( MElement *e, gmshVector & v) const;
- virtual int sizeOfR () const {return 3*_gm->getNumMeshVertices();}
- virtual int sizeOfC () const {return 3*_gm->getNumMeshVertices();}
- virtual void getComponents (int iVertex, int *components) const {
- components[0] = 3*(iVertex - 1) + 0;
- components[1] = 3*(iVertex - 1) + 1;
- components[2] = 3*(iVertex - 1) + 2;
- }
-};
-
-
-void gmshLaplaceVectorTerm:: elementMatrix ( MElement *e, gmshSmallMatrix & m) const{
- int nbNodes = e->getNumVertices();
- int integrationOrder = 2*(e->getPolynomialOrder()-1);
- int npts;
- IntPt *GP;
- double jac [3][3];
- double invjac [3][3];
- double Grads[256][3],grad[3];
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
-
- m.set_all(0.0);
-
- double FACT = _E / (1 + _nu);
- double C11 = FACT * (1 - _nu) / (1 - 2 * _nu);
- double C12 = FACT * _nu / (1 - 2 * _nu);
- double C44 = (C11 - C12) / 2;
- const double C[6][6] =
- { {C11, C12, C12, 0, 0, 0},
- {C12, C11, C12, 0, 0, 0},
- {C12, C12, C11, 0, 0, 0},
- { 0, 0, 0, C44, 0, 0},
- { 0, 0, 0, 0, C44, 0},
- { 0, 0, 0, 0, 0, C44} };
-
- gmshSmallMatrix H(6,6);
- gmshSmallMatrix B (6,3*nbNodes);
- gmshSmallMatrix BTH (3*nbNodes,6);
- gmshSmallMatrix BT (3*nbNodes,6);
- for (int i=0;i<6;i++)
- for (int j=0;j<6;j++)
- H(i,j) = C[i][j];
-
- 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);
- inv3x3 ( jac, invjac) ;
- B.set_all(0.0);
- BT.set_all(0.0);
- for (int j=0;j<nbNodes;j++){
- e->getGradShapeFunction (j,u,v,w,grad);
- Grads[j][0] = invjac[0][0] * grad[0] + invjac[0][1] * grad[1] + invjac[0][2] * grad[2];
- Grads[j][1] = invjac[1][0] * grad[0] + invjac[1][1] * grad[1] + invjac[1][2] * grad[2];
- Grads[j][2] = invjac[2][0] * grad[0] + invjac[2][1] * grad[1] + invjac[2][2] * grad[2];
- BT(j,0) = B(0,j) = Grads[j][0];
- BT(j,3) = B(3,j) = Grads[j][1];
- BT(j,4) = B(4,j) = Grads[j][2];
-
- BT(j+nbNodes,1) = B(1,j+nbNodes) = Grads[j][1];
- BT(j+nbNodes,3) = B(3,j+nbNodes) = Grads[j][0];
- BT(j+nbNodes,5) = B(5,j+nbNodes) = Grads[j][2];
-
- BT(j+2*nbNodes,2) = B(2,j+2*nbNodes) = Grads[j][2];
- BT(j+2*nbNodes,4) = B(4,j+2*nbNodes) = Grads[j][0];
- BT(j+2*nbNodes,5) = B(5,j+2*nbNodes) = Grads[j][1];
- }
- BTH.set_all(0.0);
- BTH.blas_dgemm (BT,H);
- m.blas_dgemm (BTH,B,weight*detJ,1.0);
- }
-// for (int i=0;i<12;i++){
-// for (int j=0;j<12;j++){
-// printf("%12.5E ",m(i,j));
-// }
-// printf("\n");
-// }
-// printf("\n");
-}
-
-void gmshLaplaceVectorTerm:: elementVector ( MElement *e, gmshVector & m) const{
-}
-
-/*
-Tentative to smooth the mesh according to the levelset
-*/
-
-#ifdef HAVE_GLEVELSETS
-void buildLagrangeMulipliersTerms (std::multimap<MElement*,MElement*> &cut,
- gmshSmallMatrix &C,
- gmshSmallMatrix &K){
-
- std::multimap<MElement*,MElement*>::iterator it2 = cut.begin();
- double jac[3][3];
-
- for ( ; it2 != cut.end() ; ++it2){
- MElement *e = it2->second;
- MElement *e_vol = it2->first;
-
- const double h = e->minEdge();
-
- int integrationOrder = 2*e->getPolynomialOrder();
- // Stabilization
- gmshSmallMatrix kloc(e->getNumVertices(),e->getNumVertices());
- gmshLaplaceTerm Laplace(0,.02*h*h);
- Laplace.elementMatrix (e,kloc);
-
- // Lagrange multipliers
- int nbNodes = e->getNumVertices();
- int nbNodes_vol = e_vol->getNumVertices();
- gmshSmallMatrix cloc(nbNodes,nbNodes_vol);
- int npts;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- 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);
- double xyz[3] = {0,0,0};
- double uvw[3];
- for (int k=0;k<nbNodes;k++){
- double sf;
- e->getShapeFunction (k,u,v,w,sf);
- xyz[0] += e->getVertex(k)->x() * sf;
- xyz[1] += e->getVertex(k)->y() * sf;
- xyz[2] += e->getVertex(k)->z() * sf;
- }
- e_vol->xyz2uvw(xyz,uvw);
- for (int j=0;j<nbNodes;j++){
- double sf ;e->getShapeFunction(j,u,v,w,sf);
- for (int k=0;k<nbNodes_vol;k++){
- double sf_vol;e_vol->getShapeFunction(k,uvw[0],uvw[1],uvw[2],sf_vol);
-// printf("k = %d : uvw (%g %g %g) %g %g %g %g\n",k,
-// uvw[0],uvw[1],uvw[2],
-// detJ,sf_vol,sf,weight);
- cloc(j,k) += detJ * weight * sf * sf_vol;
- }
- }
- }
-
- for (int j=0;j<nbNodes;j++){
- for (int k=0;k<nbNodes_vol;k++){
- C(e->getVertex(j)->getNum()-1,e_vol->getVertex(k)->getNum()-1)+=cloc(j,k);
- }
- for (int k=0;k<nbNodes;k++){
- K(e->getVertex(j)->getNum()-1,e->getVertex(k)->getNum()-1)+=kloc(j,k);
- }
- }
- }
-}
-#endif
-
-
-// solve a pb with lagrange multipliers
-
-int main(int argc, char **argv)
-{
- // globals are still present in Gmsh
- GmshInitialize(argc, argv);
-
- // Creation of a geometric model
- GModel *m = new GModel();
- // a Mesh is read
- m->readMSH(argv[1]);
- int nbNodes = m->getNumMeshVertices();
- int dim = m->getNumRegions() ? 3 : 2;
-
- gmshLaplaceTerm Laplace(m);
- int sP = Laplace.sizeOfR();
-
-#ifdef HAVE_GLEVELSETS
- //gLevelset *ls = new gLevelsetPlane (1,0,0,.5);
- gLevelset *ls1 = new gLevelsetSphere (-.7,.5,.5,.12);
- gLevelset *ls2 = new gLevelsetSphere (-.3,.5,.5,.12);
- // gLevelset *ls = new gLevelsetSphere (-.35,.5,.5,.2);
- std::vector<const gLevelset*> vls;
- vls.push_back(ls1);
- vls.push_back(ls2);
- gLevelsetUnion *ls = new gLevelsetUnion (vls);
- Laplace.setLevelset (ls);
- std::set<gmshTopoVertex> tVerts;
- std::multimap<MElement*,MElement*> cut;
- gmshBuildSurfaceMesh(m,ls,tVerts,cut);
- printf("%d noeuds %d triangles\n",tVerts.size(),cut.size());
- {
- FILE *fp = fopen ("cut.msh","w");
- fprintf(fp, "$MeshFormat\n2 0 8\n$EndMeshFormat\n");
- fprintf(fp, "$Nodes\n");
- fprintf(fp, "%d\n",tVerts.size());
- std::set<gmshTopoVertex>::iterator it = tVerts.begin();
- int K=1;
- for ( ; it != tVerts.end() ; ++it){
- it->getMVertex()->setNum(K);
- it->getMVertex()->setIndex (K++);
- it->getMVertex()->writeMSH(fp);
- }
- fprintf(fp, "$EndNodes\n");
- fprintf(fp, "$Elements\n");
- fprintf(fp, "%d\n",cut.size());
- std::multimap<MElement*,MElement*>::iterator it2 = cut.begin();
- for ( ; it2 != cut.end() ; ++it2){
- it2->second->writeMSH(fp,2.0);
- }
- fprintf(fp, "$EndElements\n");
- fclose(fp);
- }
- gmshSmallMatrix kstab (tVerts.size(),tVerts.size());
- gmshSmallMatrix C (tVerts.size(),m->getNumMeshVertices());
- printf("assembling lagrange terms\n");
- buildLagrangeMulipliersTerms (cut,C,kstab);
- printf("lagrange terms assembled\n");
-
- int sP1 = sP;
- sP += tVerts.size();
-#endif
- Laplace.addDirichlet (300,dim-1,0,gmshVertexEvaluator(1.0));
- Laplace.addDirichlet (200,dim-1,0,gmshVertexEvaluator(1.0));
- Laplace.addDirichlet (100,dim-1,0,gmshVertexEvaluator(1.0));
-
-#ifdef HAVE_SPARSKIT
- gmshMatrix k (sP);
-#else
- gmshMatrix k (sP,sP);
-#endif
-
- gmshVector x (sP);
- gmshVector b (sP);
-
- clock_t t1 = clock();
- Laplace.addToJacobian(k);
- Laplace.uglyDirichlet (k,b);
-
- printf("classical part assembled\n");
-
-#ifdef HAVE_GLEVELSETS
- for (int i=0;i<tVerts.size();i++)
- for (int j=0;j<tVerts.size();j++)
- addMatrix(sP1 + i,sP1 + j,kstab(i,j),k);
-
- for (int i=0;i<tVerts.size();i++){// tverts
- for (int j=0;j<sP1;j++){//totalnumv
- addMatrix(i+sP1,j,C(i,j),k);
- addMatrix(j,i+sP1,C(i,j),k);
- }
- }
-#endif
-
-
- clock_t t2 = clock();
- printf("%lf seconds for assembling the operator\n",(double)(t2-t1)/CLOCKS_PER_SEC);
- SystemSolve(k,x,b);
- clock_t t3 = clock();
- printf("%lf seconds for solving the system (%d unknowns)\n",
- (double)(t3-t2)/CLOCKS_PER_SEC,sP);
- {
- FILE *fp = fopen("res.msh","w");
- fprintf(fp, "$MeshFormat\n2 0 8\n$EndMeshFormat\n");
- fprintf(fp, "$NodeData\n");
- fprintf(fp, "1\n\"%s\"\n", "toto");
- fprintf(fp, "1\n%.16g\n", 0.0);
- fprintf(fp, "3\n%d\n%d\n%d\n", 0, 1, m->getNumMeshVertices());
- int c[100];
- for (int i=0;i<nbNodes;i++){
- Laplace.getComponents(i+1,c);
- fprintf(fp, "%d %22.15E\n",i+1, getVector(c[0],x));
- }
- fprintf(fp, "$EndNodeData\n");
- fclose (fp);
- }
-#ifdef HAVE_GLEVELSETS
- {
- FILE *fp = fopen("lagrnage.msh","w");
- fprintf(fp, "$MeshFormat\n2 0 8\n$EndMeshFormat\n");
- fprintf(fp, "$NodeData\n");
- fprintf(fp, "1\n\"%s\"\n", "toto");
- fprintf(fp, "1\n%.16g\n", 0.0);
- fprintf(fp, "3\n%d\n%d\n%d\n", 0, 1, tVerts.size());
- for (int i=0;i<tVerts.size();i++){
- fprintf(fp, "%d %22.15E\n",i+1, getVector(sP1+i,x));
- }
- fprintf(fp, "$EndNodeData\n");
- fclose (fp);
- }
-#endif
- clock_t t4 = clock();
- printf("%lf seconds for saving the solution\n",(double)(t4-t3)/CLOCKS_PER_SEC);
- GmshFinalize();
-}
-
-// solve an elasticity pb
-int main2(int argc, char **argv)
-{
- // globals are still present in Gmsh
- GmshInitialize(argc, argv);
-
- // Creation of a geometric model
- GModel *m = new GModel();
- // a Mesh is read
- m->readMSH(argv[1]);
- int nbNodes = m->getNumMeshVertices();
- int dim = m->getNumRegions() ? 3 : 2;
-
- gmshLaplaceVectorTerm Laplace(m,1.e10,.3);
- int sP = Laplace.sizeOfR();
-
- // Laplace.addDirichlet (200,dim-1,0,gmshVertexEvaluator(0.01));
- Laplace.addDirichlet (100,dim-1,0,gmshVertexEvaluator(0.0));
- Laplace.addDirichlet (100,dim-1,1,gmshVertexEvaluator(0.0));
- Laplace.addDirichlet (100,dim-1,2,gmshVertexEvaluator(0.0));
-
-
-#ifdef HAVE_SPARSKIT
- gmshMatrix k (sP);
-#else
- gmshMatrix k (sP,sP);
-#endif
- gmshVector x (sP);
- gmshVector b (sP);
-
- Laplace.addNeumann (200,dim-1,0,gmshVertexEvaluator(0.01),b);
-
- clock_t t1 = clock();
- Laplace.addToJacobian(k);
- // Laplace.addToRightHandSide(b);
- Laplace.uglyDirichlet (k,b);
- clock_t t2 = clock();
- printf("%lf seconds for assembling the operator\n",(double)(t2-t1)/CLOCKS_PER_SEC);
- SystemSolve(k,x,b);
- clock_t t3 = clock();
- // printf("%lf seconds for solving the system (%d unknowns)\n",(double)(t3-t2)/CLOCKS_PER_SEC,k.size1());
-
- FILE *fp = fopen("displacements.msh","w");
- fprintf(fp, "$MeshFormat\n2 0 8\n$EndMeshFormat\n");
- fprintf(fp, "$NodeData\n");
- fprintf(fp, "1\n\"%s\"\n", "toto");
- fprintf(fp, "1\n%.16g\n", 0.0);
- fprintf(fp, "3\n%d\n%d\n%d\n", 0, 3, m->getNumMeshVertices());
- int c[100];
- for (int i=0;i<nbNodes;i++){
- Laplace.getComponents(i+1,c);
- fprintf(fp, "%d %22.15E %22.15E %22.15E\n",i+1,
- getVector(c[0],x),getVector(c[1],x),getVector(c[2],x));
- }
- fprintf(fp, "$EndNodeData\n");
- fclose (fp);
- clock_t t4 = clock();
- printf("%lf seconds for saving the solution\n",(double)(t4-t3)/CLOCKS_PER_SEC);
- GmshFinalize();
-}
-