From af2a4e84bde774bebdcd9eb9578928828888fbda Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 3 May 2007 08:50:02 +0000
Subject: [PATCH] add test version of fourier model
---
Geo/FEdge.cpp | 25 +
Geo/FEdge.h | 41 +
Geo/FFace.cpp | 74 ++
Geo/FFace.h | 49 ++
Geo/FVertex.h | 43 +
Geo/GModel.h | 2 +-
Geo/GModelIO_F.cpp | 52 ++
Geo/Makefile | 36 +-
Parser/OpenFile.cpp | 904 +++++++++++----------
configure | 45 +-
configure.in | 30 +-
contrib/FourierModel/ContinuationPatch.cpp | 500 ++++++++++++
contrib/FourierModel/ContinuationPatch.h | 63 ++
contrib/FourierModel/Curve.cpp | 242 ++++++
contrib/FourierModel/Curve.h | 37 +
contrib/FourierModel/ExactPatch.cpp | 192 +++++
contrib/FourierModel/ExactPatch.h | 30 +
contrib/FourierModel/FM_Edge.cpp | 123 +++
contrib/FourierModel/FM_Edge.h | 37 +
contrib/FourierModel/FM_Face.cpp | 147 ++++
contrib/FourierModel/FM_Face.h | 33 +
contrib/FourierModel/FM_Info.cpp | 1 +
contrib/FourierModel/FM_Info.h | 37 +
contrib/FourierModel/FM_Reader.cpp | 108 +++
contrib/FourierModel/FM_Reader.h | 38 +
contrib/FourierModel/FM_Vertex.cpp | 2 +
contrib/FourierModel/FM_Vertex.h | 21 +
contrib/FourierModel/Main.cpp | 117 +++
contrib/FourierModel/Makefile | 55 ++
contrib/FourierModel/Message.cpp | 156 ++++
contrib/FourierModel/Message.h | 36 +
contrib/FourierModel/Patch.cpp | 84 ++
contrib/FourierModel/Patch.h | 40 +
33 files changed, 2905 insertions(+), 495 deletions(-)
create mode 100644 Geo/FEdge.cpp
create mode 100644 Geo/FEdge.h
create mode 100644 Geo/FFace.cpp
create mode 100644 Geo/FFace.h
create mode 100644 Geo/FVertex.h
create mode 100644 Geo/GModelIO_F.cpp
create mode 100644 contrib/FourierModel/ContinuationPatch.cpp
create mode 100644 contrib/FourierModel/ContinuationPatch.h
create mode 100644 contrib/FourierModel/Curve.cpp
create mode 100644 contrib/FourierModel/Curve.h
create mode 100644 contrib/FourierModel/ExactPatch.cpp
create mode 100644 contrib/FourierModel/ExactPatch.h
create mode 100644 contrib/FourierModel/FM_Edge.cpp
create mode 100644 contrib/FourierModel/FM_Edge.h
create mode 100644 contrib/FourierModel/FM_Face.cpp
create mode 100644 contrib/FourierModel/FM_Face.h
create mode 100644 contrib/FourierModel/FM_Info.cpp
create mode 100644 contrib/FourierModel/FM_Info.h
create mode 100644 contrib/FourierModel/FM_Reader.cpp
create mode 100644 contrib/FourierModel/FM_Reader.h
create mode 100644 contrib/FourierModel/FM_Vertex.cpp
create mode 100644 contrib/FourierModel/FM_Vertex.h
create mode 100644 contrib/FourierModel/Main.cpp
create mode 100644 contrib/FourierModel/Makefile
create mode 100644 contrib/FourierModel/Message.cpp
create mode 100644 contrib/FourierModel/Message.h
create mode 100644 contrib/FourierModel/Patch.cpp
create mode 100644 contrib/FourierModel/Patch.h
diff --git a/Geo/FEdge.cpp b/Geo/FEdge.cpp
new file mode 100644
index 0000000000..545c57024d
--- /dev/null
+++ b/Geo/FEdge.cpp
@@ -0,0 +1,25 @@
+#include "FEdge.h"
+
+#if defined(HAVE_FOURIER_MODEL)
+
+Range<double> FEdge::parBounds(int i) const
+{
+ return(Range<double>(0.,1.));
+}
+
+GPoint FEdge::point(double p) const
+{
+ double x, y, z;
+
+ edge->F(p,x,y,z);
+ return GPoint(x,y,z);
+}
+
+double FEdge::parFromPoint(const SPoint3 &pt) const
+{
+ double t;
+ edge->Inverse(pt.x(),pt.y(),pt.z(),t);
+ return t;
+}
+
+#endif
diff --git a/Geo/FEdge.h b/Geo/FEdge.h
new file mode 100644
index 0000000000..ba59a9378c
--- /dev/null
+++ b/Geo/FEdge.h
@@ -0,0 +1,41 @@
+#ifndef _F_EDGE_H_
+#define _F_EDGE_H_
+
+#include "GEdge.h"
+#include "GModel.h"
+#include "FVertex.h"
+#include "Range.h"
+
+#if defined(HAVE_FOURIER_MODEL)
+
+#include "FM_Edge.h"
+
+class FEdge : public GEdge {
+ protected:
+ FM_Edge* edge;
+ public:
+ FEdge(GModel *model, FM_Edge* edge_, int tag, GVertex *v1, GVertex *v2) :
+ GEdge(model, tag, v1, v2), edge(edge_) {}
+ virtual ~FEdge() {}
+ double period() const { throw ; }
+ virtual bool periodic(int dim=0) const { return false; }
+ virtual Range<double> parBounds(int i) const;
+ virtual GeomType geomType() const { throw; }
+ virtual bool degenerate(int) const { return false; }
+ virtual bool continuous(int dim) const { return true; }
+ virtual GPoint point(double p) const;
+ virtual GPoint closestPoint(const SPoint3 & queryPoint) { throw; }
+ virtual int containsPoint(const SPoint3 &pt) const { throw; }
+ virtual int containsParam(double pt) const { throw; }
+ virtual SVector3 firstDer(double par) const { throw; }
+ virtual SPoint2 reparamOnFace(GFace *face, double epar, int dir) const
+ { throw; }
+ virtual double parFromPoint(const SPoint3 &pt) const;
+ virtual int minimumMeshSegments () const { throw; }
+ virtual int minimumDrawSegments () const { throw; }
+ ModelType getNativeType() const { return FourierModel; }
+};
+
+#endif
+
+#endif
diff --git a/Geo/FFace.cpp b/Geo/FFace.cpp
new file mode 100644
index 0000000000..6640975444
--- /dev/null
+++ b/Geo/FFace.cpp
@@ -0,0 +1,74 @@
+#include "FVertex.h"
+#include "FFace.h"
+
+#if defined(HAVE_FOURIER_MODEL)
+
+FFace::FFace(GModel *m, FM_Face *face_, int tag) : GFace(m,tag), face(face_)
+{
+ int curCorner = 0;
+ for (int i=0;i<face->GetNumEdges();i++) {
+ GVertex* p0 = new FVertex(m,curCorner,face->GetEdge(i)->GetStartPoint());
+ curCorner = (curCorner + 1) % face->GetNumEdges();
+ GVertex* p1 = new FVertex(m,curCorner,face->GetEdge(i)->GetEndPoint());
+ l_edges.push_back(new FEdge(m,face->GetEdge(i),i,p0,p1));
+ l_dirs.push_back(1);
+ }
+}
+
+Range<double> FFace::parBounds(int i) const
+{
+ return Range<double>(0.,1.);
+}
+
+GPoint FFace::point(double par1, double par2) const
+{
+ double pp[2] = {par1,par2};
+ double x,y,z;
+ face->F(par1,par2,x,y,z);
+ return GPoint(x, y, z, this, pp);
+}
+
+GPoint FFace::point(const SPoint2 &pt) const
+{
+ return point(pt[0], pt[1]);
+}
+
+SPoint2 FFace::parFromPoint(const SPoint3 &p) const
+{
+ double u, v, x, y, z;
+ x = p.x(); y = p.y(); z = p.z();
+ face->Inverse(x,y,z,u,v);
+ return SPoint2(u, v);
+}
+
+GPoint FFace::closestPoint(const SPoint3 & queryPoint) const
+{
+ throw;
+}
+
+int FFace::containsPoint(const SPoint3 &pt) const
+{
+ throw;
+}
+
+int FFace::containsParam(const SPoint2 &pt) const
+{
+ const double tol = 1.e-6;
+ if(pt[0] < 0. - tol || pt[0] > 1. + tol) return 0;
+ if(pt[1] < 0. - tol || pt[1] > 1. + tol) return 0;
+ return 1;
+}
+
+SVector3 FFace::normal(const SPoint2 ¶m) const
+{
+ double x,y,z;
+ face->GetUnitNormal(param[0],param[1],x,y,z);
+ return SVector3(x, y, z);
+}
+
+GEntity::GeomType FFace::geomType() const
+{
+ return GEntity::ParametricSurface;
+}
+
+#endif
diff --git a/Geo/FFace.h b/Geo/FFace.h
new file mode 100644
index 0000000000..5c6bb11b8f
--- /dev/null
+++ b/Geo/FFace.h
@@ -0,0 +1,49 @@
+#ifndef _F_FACE_H_
+#define _F_FACE_H_
+
+#include "GFace.h"
+#include "GModel.h"
+#include "Range.h"
+#include "FEdge.h"
+
+#if defined(HAVE_FOURIER_MODEL)
+
+#include "FM_Face.h"
+
+class FFace : public GFace {
+ protected:
+ FM_Face *face;
+ bool _periodic[2];
+ public:
+ FFace(GModel *m, FM_Face *face_, int tag);
+
+ virtual ~FFace() {}
+
+ Range<double> parBounds(int i) const;
+ virtual int paramDegeneracies(int dir, double *par) { return 0; }
+ virtual GPoint point(double par1, double par2) const;
+ virtual GPoint point(const SPoint2 &pt) const;
+ virtual SPoint2 parFromPoint(const SPoint3 &p) const;
+ virtual GPoint closestPoint(const SPoint3 & queryPoint) const;
+ virtual int containsPoint(const SPoint3 &pt) const;
+ virtual int containsParam(const SPoint2 &pt) const;
+ virtual SVector3 normal(const SPoint2 ¶m) const;
+ virtual GEntity::GeomType geomType() const;
+ virtual Pair<SVector3,SVector3>
+ firstDer(const SPoint2 ¶m) const {throw;}
+ virtual double * nthDerivative
+ (const SPoint2 ¶m, int n, double *array) const {throw;}
+ virtual int geomDirection() const { return 1; }
+ virtual bool continuous(int dim) const { return true; }
+ virtual bool periodic(int dim) const { return false; }
+ virtual bool degenerate(int dim) const { return false; }
+ virtual double period(int dir) const {throw;}
+ ModelType getNativeType() const { return FourierModel; }
+ void * getNativePtr() const {throw;}
+ virtual bool surfPeriodic(int dim) const {throw;}
+};
+
+#endif
+
+#endif
+
diff --git a/Geo/FVertex.h b/Geo/FVertex.h
new file mode 100644
index 0000000000..481e16fe03
--- /dev/null
+++ b/Geo/FVertex.h
@@ -0,0 +1,43 @@
+#ifndef _F_VERTEX_H_
+#define _F_VERTEX_H_
+
+#include "GModel.h"
+#include "GVertex.h"
+
+#if defined(HAVE_FOURIER_MODEL)
+
+#include "FM_Vertex.h"
+
+class FVertex : public GVertex {
+ protected:
+ FM_Vertex* v;
+
+ public:
+ FVertex(GModel *m, int num, FM_Vertex* _v) : GVertex(m, num), v(_v)
+ {
+ mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
+ }
+ virtual ~FVertex() {}
+ virtual GPoint point() const
+ {
+ return GPoint(x(),y(),z());
+ }
+ virtual double x() const
+ {
+ return v->GetX();
+ }
+ virtual double y() const
+ {
+ return v->GetY();
+ }
+ virtual double z() const
+ {
+ return v->GetZ();
+ }
+ ModelType getNativeType() const { return FourierModel; }
+ virtual SPoint2 reparamOnFace ( GFace *gf , int) const { throw; }
+};
+
+#endif
+
+#endif
diff --git a/Geo/GModel.h b/Geo/GModel.h
index e03636bb6d..20487dd475 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -167,7 +167,7 @@ class GModel
int writeGEO(const std::string &name, bool printLabels=true);
// Fourier model
- int readFourier(const std::string &name);
+ int readF(const std::string &name);
// OCC model
int readOCCBREP(const std::string &name);
diff --git a/Geo/GModelIO_F.cpp b/Geo/GModelIO_F.cpp
new file mode 100644
index 0000000000..c096b4a63f
--- /dev/null
+++ b/Geo/GModelIO_F.cpp
@@ -0,0 +1,52 @@
+#include <string>
+#include "GModel.h"
+#include "Message.h"
+#include "Context.h"
+#include "Views.h"
+#include "FFace.h"
+#include "meshGFace.h"
+
+#if defined(HAVE_FOURIER_MODEL)
+
+#include "FM_Reader.h"
+
+extern Context_T CTX;
+
+class meshGmsh{
+public:
+ void operator() (GFace *gf)
+ {
+ FFace *ff = dynamic_cast<FFace*>(gf);
+ if(!ff) {
+ Msg(GERROR, "face %d is not Fourier", gf->tag());
+ return;
+ }
+ meshGFace mesh;
+ mesh(ff);
+ }
+};
+
+int GModel::readF(const std::string &fn)
+{
+ CTX.terminal = 1;
+
+ FM_Reader reader(fn.c_str());
+
+ for (int i=0;i<reader.GetNumFaces();i++) {
+ add(new FFace(this, reader.GetFace(i), i));
+ }
+ Msg(INFO, "Fourier model created: %d patches", reader.GetNumFaces());
+
+ return 1;
+}
+
+#else
+
+int GModel::readF(const std::string &fn)
+{
+ Msg(GERROR,"Gmsh has to be compiled with Fourier Model support to load %s",
+ fn.c_str());
+ return 0;
+}
+
+#endif
diff --git a/Geo/Makefile b/Geo/Makefile
index a4ea4b0aac..d818f5efe8 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.143 2007-04-23 08:04:16 geuzaine Exp $
+# $Id: Makefile,v 1.144 2007-05-03 08:50:01 geuzaine Exp $
#
# Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
#
@@ -30,12 +30,13 @@ SRC = GEntity.cpp\
GVertex.cpp GEdge.cpp GEdgeLoop.cpp GFace.cpp GRegion.cpp\
gmshVertex.cpp gmshEdge.cpp gmshFace.cpp gmshRegion.cpp gmshSurface.cpp\
OCCVertex.cpp OCCEdge.cpp OCCFace.cpp OCCRegion.cpp\
+ FEdge.cpp FFace.cpp\
projectionFace.cpp\
GModel.cpp\
GModelIO_Geo.cpp\
GModelIO_Mesh.cpp\
- GModelIO_Fourier.cpp\
GModelIO_OCC.cpp\
+ GModelIO_F.cpp\
GModelIO_CGNS.cpp\
GModelIO_MED.cpp\
ExtrudeParams.cpp \
@@ -160,6 +161,18 @@ OCCRegion.o: OCCRegion.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
ExtrudeParams.h ../Common/SmoothData.h GFace.h GEdgeLoop.h Pair.h \
GRegion.h OCCVertex.h OCCIncludes.h OCCEdge.h OCCFace.h OCCRegion.h \
../Common/Message.h
+FEdge.o: FEdge.cpp FEdge.h GEdge.h GEntity.h Range.h SPoint3.h \
+ SBoundingBox3d.h ../Common/GmshDefines.h GVertex.h MVertex.h GPoint.h \
+ SPoint2.h SVector3.h MElement.h MEdge.h ../Common/Hash.h MFace.h \
+ ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
+ ExtrudeParams.h ../Common/SmoothData.h GModel.h GFace.h GEdgeLoop.h \
+ Pair.h GRegion.h FVertex.h
+FFace.o: FFace.cpp FVertex.h GModel.h GVertex.h GEntity.h Range.h \
+ SPoint3.h SBoundingBox3d.h ../Common/GmshDefines.h MVertex.h GPoint.h \
+ SPoint2.h GEdge.h SVector3.h MElement.h MEdge.h ../Common/Hash.h \
+ MFace.h ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
+ ExtrudeParams.h ../Common/SmoothData.h GFace.h GEdgeLoop.h Pair.h \
+ GRegion.h FFace.h FEdge.h
projectionFace.o: projectionFace.cpp projectionFace.h GFace.h GPoint.h \
GEntity.h Range.h SPoint3.h SBoundingBox3d.h ../Common/GmshDefines.h \
GEdgeLoop.h GEdge.h GVertex.h MVertex.h SPoint2.h SVector3.h MElement.h \
@@ -195,16 +208,6 @@ GModelIO_Mesh.o: GModelIO_Mesh.cpp ../Common/Message.h \
../Common/SmoothData.h GFace.h GEdgeLoop.h Pair.h GRegion.h \
gmshRegion.h Geo.h gmshSurface.h ../DataStr/Tree.h ../DataStr/avl.h \
gmshFace.h gmshVertex.h gmshEdge.h
-GModelIO_Fourier.o: GModelIO_Fourier.cpp GModel.h GVertex.h GEntity.h \
- Range.h SPoint3.h SBoundingBox3d.h ../Common/GmshDefines.h MVertex.h \
- GPoint.h SPoint2.h GEdge.h SVector3.h MElement.h MEdge.h \
- ../Common/Hash.h MFace.h ../Numeric/Numeric.h ../Common/Context.h \
- ../DataStr/List.h ExtrudeParams.h ../Common/SmoothData.h GFace.h \
- GEdgeLoop.h Pair.h GRegion.h fourierFace.h gmshFace.h Geo.h \
- gmshSurface.h ../DataStr/Tree.h ../DataStr/avl.h gmshVertex.h \
- ../Common/Message.h ../Common/Views.h ../Common/ColorTable.h \
- ../Common/VertexArray.h ../Common/SmoothData.h \
- ../Common/AdaptiveViews.h ../Common/GmshMatrix.h
GModelIO_OCC.o: GModelIO_OCC.cpp GModelIO_OCC.h GModel.h GVertex.h \
GEntity.h Range.h SPoint3.h SBoundingBox3d.h ../Common/GmshDefines.h \
MVertex.h GPoint.h SPoint2.h GEdge.h SVector3.h MElement.h MEdge.h \
@@ -212,6 +215,15 @@ GModelIO_OCC.o: GModelIO_OCC.cpp GModelIO_OCC.h GModel.h GVertex.h \
../DataStr/List.h ExtrudeParams.h ../Common/SmoothData.h GFace.h \
GEdgeLoop.h Pair.h GRegion.h OCCIncludes.h ../Common/Message.h \
OCCVertex.h OCCEdge.h OCCFace.h OCCRegion.h
+GModelIO_F.o: GModelIO_F.cpp GModel.h GVertex.h GEntity.h Range.h \
+ SPoint3.h SBoundingBox3d.h ../Common/GmshDefines.h MVertex.h GPoint.h \
+ SPoint2.h GEdge.h SVector3.h MElement.h MEdge.h ../Common/Hash.h \
+ MFace.h ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
+ ExtrudeParams.h ../Common/SmoothData.h GFace.h GEdgeLoop.h Pair.h \
+ GRegion.h ../Common/Message.h ../Common/Views.h ../Common/ColorTable.h \
+ ../Common/VertexArray.h ../Common/SmoothData.h \
+ ../Common/AdaptiveViews.h ../Common/GmshMatrix.h FFace.h FEdge.h \
+ FVertex.h ../Mesh/meshGFace.h
GModelIO_CGNS.o: GModelIO_CGNS.cpp GModel.h GVertex.h GEntity.h Range.h \
SPoint3.h SBoundingBox3d.h ../Common/GmshDefines.h MVertex.h GPoint.h \
SPoint2.h GEdge.h SVector3.h MElement.h MEdge.h ../Common/Hash.h \
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index 7474acff9d..0024d2836b 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,450 +1,454 @@
-// $Id: OpenFile.cpp,v 1.145 2007-03-16 10:03:40 remacle Exp $
-//
-// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
-//
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 2 of the License, or
-// (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
-// USA.
-//
-// Please report all bugs and problems to <gmsh@geuz.org>.
-
-#if defined(__CYGWIN__)
-#include <sys/cygwin.h>
-#endif
-
-#include "Gmsh.h"
-#include "Geo.h"
-#include "GModel.h"
-#include "Numeric.h"
-#include "Context.h"
-#include "Parser.h"
-#include "OpenFile.h"
-#include "CommandLine.h"
-#include "Views.h"
-#include "ReadImg.h"
-#include "OS.h"
-#include "HighOrder.h"
-
-#if defined(HAVE_FLTK)
-#include "GmshUI.h"
-#include "Draw.h"
-#include "SelectBuffer.h"
-#include "GUI.h"
-extern GUI *WID;
-void UpdateViewsInGUI();
-#endif
-
-extern Mesh *THEM;
-extern GModel *GMODEL;
-extern Context_T CTX;
-
-void FixRelativePath(char *in, char *out){
- if(in[0] == '/' || in[0] == '\\' || (strlen(in)>2 && in[1] == ':')){
- // do nothing: 'in' is an absolute path
- strcpy(out, in);
- }
- else{
- // append 'in' to the path of the parent file
- strcpy(out, yyname);
- int i = strlen(out)-1 ;
- while(i >= 0 && yyname[i] != '/' && yyname[i] != '\\') i-- ;
- out[i+1] = '\0';
- strcat(out, in);
- }
-}
-
-void FixWindowsPath(char *in, char *out){
-#if defined(__CYGWIN__)
- cygwin_conv_to_win32_path(in, out);
-#else
- strcpy(out, in);
-#endif
-}
-
-void SplitFileName(char *name, char *base, char *ext)
-{
- strcpy(base, name);
- strcpy(ext, "");
- for(int i = strlen(name)-1; i >= 0; i--){
- if(name[i] == '.'){
- strcpy(ext, &name[i]);
- base[i] = '\0';
- break;
- }
- }
-}
-
-static void FinishUpBoundingBox()
-{
- double range[3];
-
- for(int i = 0; i < 3; i++){
- CTX.cg[i] = 0.5 * (CTX.min[i] + CTX.max[i]);
- range[i] = CTX.max[i] - CTX.min[i];
- }
-
- if(range[0] == 0. && range[1] == 0. && range[2] == 0.) {
- CTX.min[0] -= 1.; CTX.min[1] -= 1.;
- CTX.max[0] += 1.; CTX.max[1] += 1.;
- CTX.lc = 1.;
- }
- else if(range[0] == 0. && range[1] == 0.) {
- CTX.lc = range[2];
- CTX.min[0] -= CTX.lc; CTX.min[1] -= CTX.lc;
- CTX.max[0] += CTX.lc; CTX.max[1] += CTX.lc;
- }
- else if(range[0] == 0. && range[2] == 0.) {
- CTX.lc = range[1];
- CTX.min[0] -= CTX.lc; CTX.max[0] += CTX.lc;
- }
- else if(range[1] == 0. && range[2] == 0.) {
- CTX.lc = range[0];
- CTX.min[1] -= CTX.lc; CTX.max[1] += CTX.lc;
- }
- else if(range[0] == 0.) {
- CTX.lc = sqrt(DSQR(range[1]) + DSQR(range[2]));
- CTX.min[0] -= CTX.lc; CTX.max[0] += CTX.lc;
- }
- else if(range[1] == 0.) {
- CTX.lc = sqrt(DSQR(range[0]) + DSQR(range[2]));
- CTX.min[1] -= CTX.lc; CTX.max[1] += CTX.lc;
- }
- else if(range[2] == 0.) {
- CTX.lc = sqrt(DSQR(range[0]) + DSQR(range[1]));
- }
- else {
- CTX.lc = sqrt(DSQR(range[0]) + DSQR(range[1]) + DSQR(range[2]));
- }
-}
-
-void SetBoundingBox(double xmin, double xmax,
- double ymin, double ymax,
- double zmin, double zmax)
-{
- CTX.min[0] = xmin; CTX.max[0] = xmax;
- CTX.min[1] = ymin; CTX.max[1] = ymax;
- CTX.min[2] = zmin; CTX.max[2] = zmax;
- FinishUpBoundingBox();
-}
-
-void SetBoundingBox(void)
-{
- if(CTX.forced_bbox) return;
-
- SBoundingBox3d bb;
-
- bb = GMODEL->bounds();
-
- if(bb.empty() && List_Nbr(CTX.post.list)) {
- for(int i = 0; i < List_Nbr(CTX.post.list); i++){
- Post_View *v = *(Post_View **)List_Pointer(CTX.post.list, i);
- if(fabs(v->BBox[0]) != VAL_INF &&
- fabs(v->BBox[2]) != VAL_INF &&
- fabs(v->BBox[4]) != VAL_INF)
- bb += SPoint3(v->BBox[0], v->BBox[2], v->BBox[4]);
- if(fabs(v->BBox[1]) != VAL_INF &&
- fabs(v->BBox[3]) != VAL_INF &&
- fabs(v->BBox[5]) != VAL_INF)
- bb += SPoint3(v->BBox[1], v->BBox[3], v->BBox[5]);
- }
- }
-
- if(bb.empty()){
- bb += SPoint3(-1., -1., -1.);
- bb += SPoint3(1., 1., 1.);
- }
-
- CTX.min[0] = bb.min().x(); CTX.max[0] = bb.max().x();
- CTX.min[1] = bb.min().y(); CTX.max[1] = bb.max().y();
- CTX.min[2] = bb.min().z(); CTX.max[2] = bb.max().z();
- FinishUpBoundingBox();
-}
-
-// FIXME: this is necessary for now to have an approximate CTX.lc
-// *while* parsing input files (it's important since some of the
-// geometrical operations use a tolerance that depends on
-// CTX.lc). This will be removed once the new database is filled
-// directly during the parsing step
-static SBoundingBox3d temp_bb;
-
-void ResetTemporaryBoundingBox()
-{
- temp_bb.reset();
-}
-
-void AddToTemporaryBoundingBox(double x, double y, double z)
-{
- temp_bb += SPoint3(x, y, z);
- CTX.min[0] = temp_bb.min().x(); CTX.max[0] = temp_bb.max().x();
- CTX.min[1] = temp_bb.min().y(); CTX.max[1] = temp_bb.max().y();
- CTX.min[2] = temp_bb.min().z(); CTX.max[2] = temp_bb.max().z();
- FinishUpBoundingBox();
-}
-
-int ParseFile(char *f, int close, int warn_if_missing)
-{
- char yyname_old[256], tmp[256];
- FILE *yyin_old, *fp;
- int yylineno_old, yyerrorstate_old, numviews_old;
-
- // add 'b' for pure Windows programs: opening in text mode messes up
- // fsetpos/fgetpos (used e.g. for user-defined functions)
- if(!(fp = fopen(f, "rb"))){
- if(warn_if_missing) Msg(WARNING, "Unable to open file '%s'", f);
- return 0;
- }
-
- strncpy(yyname_old, yyname, 255);
- yyin_old = yyin;
- yyerrorstate_old = yyerrorstate;
- yylineno_old = yylineno;
- numviews_old = List_Nbr(CTX.post.list);
-
- strncpy(yyname, f, 255);
- yyin = fp;
- yyerrorstate = 0;
- yylineno = 1;
-
- fpos_t position;
- fgetpos(yyin, &position);
- fgets(tmp, sizeof(tmp), yyin);
- fsetpos(yyin, &position);
-
- while(!feof(yyin)){
- yyparse();
- if(yyerrorstate > 20){
- Msg(GERROR, "Too many errors: aborting...");
- force_yyflush();
- break;
- }
- }
-
- if(close)
- fclose(yyin);
-
- strncpy(yyname, yyname_old, 255);
- yyin = yyin_old;
- yyerrorstate = yyerrorstate_old;
- yylineno = yylineno_old;
-
- if(List_Nbr(CTX.post.list) != numviews_old){
-#if defined(HAVE_FLTK)
- UpdateViewsInGUI();
-#endif
- }
-
- return 1;
-}
-
-void ParseString(char *str)
-{
- if(!str) return;
- FILE *fp;
- if((fp = fopen(CTX.tmp_filename_fullpath, "w"))) {
- fprintf(fp, str);
- fprintf(fp, "\n");
- fclose(fp);
- ParseFile(CTX.tmp_filename_fullpath, 1);
- if(GMODEL) GMODEL->importTHEM();
- }
-}
-
-void SetProjectName(char *name)
-{
- char base[356], ext[256];
- SplitFileName(name, base, ext);
-
- strncpy(CTX.filename, name, 255);
- strncpy(CTX.base_filename, base, 255);
-
-#if defined(HAVE_FLTK)
- if(!CTX.batch) WID->set_title(CTX.filename);
-#endif
-}
-
-int MergeFile(char *name, int warn_if_missing)
-{
-#if defined(HAVE_FOURIER_MODEL)
- if(!strcmp(name, "falcon") || !strcmp(name, "ship")){
- GMODEL->readFourier(name);
- SetBoundingBox();
- CTX.mesh.changed = ENT_ALL;
- return 1;
- }
-#endif
-
- // added 'b' for pure Windows programs, since some of these files
- // contain binary data
- FILE *fp = fopen(name, "rb");
- if(!fp){
- if(warn_if_missing) Msg(WARNING, "Unable to open file '%s'", name);
- return 0;
- }
-
- char header[256];
- fgets(header, sizeof(header), fp);
- fclose(fp);
-
- Msg(STATUS2, "Reading '%s'", name);
-
- char ext[256], base[256];
- SplitFileName(name, base, ext);
-
-#if defined(HAVE_FLTK)
- if(!CTX.batch) {
- if(!strcmp(ext, ".gz")) {
- // the real solution would be to rewrite all our I/O functions in
- // terms of gzFile, but until then, this is better than nothing
- if(fl_choice("File '%s' is in gzip format.\n\nDo you want to uncompress it?",
- "Cancel", "Uncompress", NULL, name)){
- char tmp[256];
- sprintf(tmp, "gunzip -c %s > %s", name, base);
- SystemCall(tmp);
- if(!strcmp(CTX.filename, name)) // this is the project file
- SetProjectName(base);
- return MergeFile(base);
- }
- }
- }
-#endif
-
- CTX.geom.draw = 0; // don't try to draw the model while reading
- int status = 0;
- if(!strcmp(ext, ".stl") || !strcmp(ext, ".STL")){
- status = GMODEL->readSTL(name, CTX.geom.tolerance);
- }
- else if(!strcmp(ext, ".brep") || !strcmp(ext, ".rle") ||
- !strcmp(ext, ".brp") || !strcmp(ext, ".BRP")){
- GMODEL->readOCCBREP(std::string(name));
- }
- else if(!strcmp(ext, ".iges") || !strcmp(ext, ".IGES") ||
- !strcmp(ext, ".igs") || !strcmp(ext, ".IGS")){
- GMODEL->readOCCIGES(std::string(name));
- }
- else if(!strcmp(ext, ".step") || !strcmp(ext, ".STEP") ||
- !strcmp(ext, ".stp") || !strcmp(ext, ".STP")){
- GMODEL->readOCCSTEP(std::string(name));
- }
- else if(!strcmp(ext, ".unv") || !strcmp(ext, ".UNV")){
- status = GMODEL->readUNV(name);
- }
- else if(!strcmp(ext, ".wrl") || !strcmp(ext, ".WRL") ||
- !strcmp(ext, ".vrml") || !strcmp(ext, ".VRML") ||
- !strcmp(ext, ".iv") || !strcmp(ext, ".IV")){
- status = GMODEL->readVRML(name);
- }
- else if(!strcmp(ext, ".mesh") || !strcmp(ext, ".MESH")){
- status = GMODEL->readMESH(name);
- }
- else if(!strcmp(ext, ".bdf") || !strcmp(ext, ".BDF") ||
- !strcmp(ext, ".nas") || !strcmp(ext, ".NAS")){
- status = GMODEL->readBDF(name);
- }
- else if(!strcmp(ext, ".p3d") || !strcmp(ext, ".P3D")){
- status = GMODEL->readP3D(name);
- }
-#if defined(HAVE_FLTK)
- else if(!strcmp(ext, ".pnm") || !strcmp(ext, ".PNM") ||
- !strcmp(ext, ".pbm") || !strcmp(ext, ".PBM") ||
- !strcmp(ext, ".pgm") || !strcmp(ext, ".PGM") ||
- !strcmp(ext, ".ppm") || !strcmp(ext, ".PPM")) {
- status = read_pnm(name);
- }
- else if(!strcmp(ext, ".bmp") || !strcmp(ext, ".BMP")) {
- status = read_bmp(name);
- }
-#if defined(HAVE_LIBJPEG)
- else if(!strcmp(ext, ".jpg") || !strcmp(ext, ".JPG") ||
- !strcmp(ext, ".jpeg") || !strcmp(ext, ".JPEG")) {
- status = read_jpeg(name);
- }
-#endif
-#if defined(HAVE_LIBPNG)
- else if(!strcmp(ext, ".png") || !strcmp(ext, ".PNG")) {
- status = read_png(name);
- }
-#endif
-#endif
- else {
- CTX.geom.draw = 1;
- if(!strncmp(header, "$PTS", 4) || !strncmp(header, "$NO", 3) ||
- !strncmp(header, "$PARA", 5) || !strncmp(header, "$ELM", 4) ||
- !strncmp(header, "$MeshFormat", 11)) {
- status = GMODEL->readMSH(name);
- }
- else if(!strncmp(header, "$PostFormat", 11) ||
- !strncmp(header, "$View", 5)) {
- status = ReadView(name);
- }
- else {
- status = GMODEL->readGEO(name);
- }
- }
-
- SetBoundingBox();
-
- CTX.geom.draw = 1;
- CTX.mesh.changed = ENT_ALL;
-
- checkHighOrderTriangles ( GMODEL );
-
- Msg(STATUS2, "Read '%s'", name);
- return status;
-}
-
-void OpenProject(char *name)
-{
- if(CTX.threads_lock) {
- Msg(INFO, "I'm busy! Ask me that later...");
- return;
- }
- CTX.threads_lock = 1;
-
- GMODEL->destroy();
- THEM->destroy();
-
- // Initialize pseudo random mesh generator to the same seed
- srand(1);
-
- // temporary hack until we fill GMODEL on the fly during parsing
- ResetTemporaryBoundingBox();
-
- SetProjectName(name);
- MergeFile(name);
-
- CTX.threads_lock = 0;
-
-#if defined(HAVE_FLTK)
- if(!CTX.batch) WID->reset_visibility();
- ZeroHighlight();
-#endif
-}
-
-void OpenProjectMacFinder(const char *filename)
-{
- static int first = 1;
- if(first || CTX.batch){
- // just copy the filename: it will be opened when the GUI is ready
- // in main()
- strncpy(CTX.filename, filename, 255);
- first = 0;
- }
- else{
- OpenProject((char*)filename);
-#if defined(HAVE_FLTK)
- Draw();
-#endif
- }
-}
-
+// $Id: OpenFile.cpp,v 1.146 2007-05-03 08:50:02 geuzaine Exp $
+//
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#if defined(__CYGWIN__)
+#include <sys/cygwin.h>
+#endif
+
+#include "Gmsh.h"
+#include "Geo.h"
+#include "GModel.h"
+#include "Numeric.h"
+#include "Context.h"
+#include "Parser.h"
+#include "OpenFile.h"
+#include "CommandLine.h"
+#include "Views.h"
+#include "ReadImg.h"
+#include "OS.h"
+#include "HighOrder.h"
+
+#if defined(HAVE_FLTK)
+#include "GmshUI.h"
+#include "Draw.h"
+#include "SelectBuffer.h"
+#include "GUI.h"
+extern GUI *WID;
+void UpdateViewsInGUI();
+#endif
+
+extern Mesh *THEM;
+extern GModel *GMODEL;
+extern Context_T CTX;
+
+void FixRelativePath(char *in, char *out){
+ if(in[0] == '/' || in[0] == '\\' || (strlen(in)>2 && in[1] == ':')){
+ // do nothing: 'in' is an absolute path
+ strcpy(out, in);
+ }
+ else{
+ // append 'in' to the path of the parent file
+ strcpy(out, yyname);
+ int i = strlen(out)-1 ;
+ while(i >= 0 && yyname[i] != '/' && yyname[i] != '\\') i-- ;
+ out[i+1] = '\0';
+ strcat(out, in);
+ }
+}
+
+void FixWindowsPath(char *in, char *out){
+#if defined(__CYGWIN__)
+ cygwin_conv_to_win32_path(in, out);
+#else
+ strcpy(out, in);
+#endif
+}
+
+void SplitFileName(char *name, char *base, char *ext)
+{
+ strcpy(base, name);
+ strcpy(ext, "");
+ for(int i = strlen(name)-1; i >= 0; i--){
+ if(name[i] == '.'){
+ strcpy(ext, &name[i]);
+ base[i] = '\0';
+ break;
+ }
+ }
+}
+
+static void FinishUpBoundingBox()
+{
+ double range[3];
+
+ for(int i = 0; i < 3; i++){
+ CTX.cg[i] = 0.5 * (CTX.min[i] + CTX.max[i]);
+ range[i] = CTX.max[i] - CTX.min[i];
+ }
+
+ if(range[0] == 0. && range[1] == 0. && range[2] == 0.) {
+ CTX.min[0] -= 1.; CTX.min[1] -= 1.;
+ CTX.max[0] += 1.; CTX.max[1] += 1.;
+ CTX.lc = 1.;
+ }
+ else if(range[0] == 0. && range[1] == 0.) {
+ CTX.lc = range[2];
+ CTX.min[0] -= CTX.lc; CTX.min[1] -= CTX.lc;
+ CTX.max[0] += CTX.lc; CTX.max[1] += CTX.lc;
+ }
+ else if(range[0] == 0. && range[2] == 0.) {
+ CTX.lc = range[1];
+ CTX.min[0] -= CTX.lc; CTX.max[0] += CTX.lc;
+ }
+ else if(range[1] == 0. && range[2] == 0.) {
+ CTX.lc = range[0];
+ CTX.min[1] -= CTX.lc; CTX.max[1] += CTX.lc;
+ }
+ else if(range[0] == 0.) {
+ CTX.lc = sqrt(DSQR(range[1]) + DSQR(range[2]));
+ CTX.min[0] -= CTX.lc; CTX.max[0] += CTX.lc;
+ }
+ else if(range[1] == 0.) {
+ CTX.lc = sqrt(DSQR(range[0]) + DSQR(range[2]));
+ CTX.min[1] -= CTX.lc; CTX.max[1] += CTX.lc;
+ }
+ else if(range[2] == 0.) {
+ CTX.lc = sqrt(DSQR(range[0]) + DSQR(range[1]));
+ }
+ else {
+ CTX.lc = sqrt(DSQR(range[0]) + DSQR(range[1]) + DSQR(range[2]));
+ }
+}
+
+void SetBoundingBox(double xmin, double xmax,
+ double ymin, double ymax,
+ double zmin, double zmax)
+{
+ CTX.min[0] = xmin; CTX.max[0] = xmax;
+ CTX.min[1] = ymin; CTX.max[1] = ymax;
+ CTX.min[2] = zmin; CTX.max[2] = zmax;
+ FinishUpBoundingBox();
+}
+
+void SetBoundingBox(void)
+{
+ if(CTX.forced_bbox) return;
+
+ SBoundingBox3d bb;
+
+ bb = GMODEL->bounds();
+
+ if(bb.empty() && List_Nbr(CTX.post.list)) {
+ for(int i = 0; i < List_Nbr(CTX.post.list); i++){
+ Post_View *v = *(Post_View **)List_Pointer(CTX.post.list, i);
+ if(fabs(v->BBox[0]) != VAL_INF &&
+ fabs(v->BBox[2]) != VAL_INF &&
+ fabs(v->BBox[4]) != VAL_INF)
+ bb += SPoint3(v->BBox[0], v->BBox[2], v->BBox[4]);
+ if(fabs(v->BBox[1]) != VAL_INF &&
+ fabs(v->BBox[3]) != VAL_INF &&
+ fabs(v->BBox[5]) != VAL_INF)
+ bb += SPoint3(v->BBox[1], v->BBox[3], v->BBox[5]);
+ }
+ }
+
+ if(bb.empty()){
+ bb += SPoint3(-1., -1., -1.);
+ bb += SPoint3(1., 1., 1.);
+ }
+
+ CTX.min[0] = bb.min().x(); CTX.max[0] = bb.max().x();
+ CTX.min[1] = bb.min().y(); CTX.max[1] = bb.max().y();
+ CTX.min[2] = bb.min().z(); CTX.max[2] = bb.max().z();
+ FinishUpBoundingBox();
+}
+
+// FIXME: this is necessary for now to have an approximate CTX.lc
+// *while* parsing input files (it's important since some of the
+// geometrical operations use a tolerance that depends on
+// CTX.lc). This will be removed once the new database is filled
+// directly during the parsing step
+static SBoundingBox3d temp_bb;
+
+void ResetTemporaryBoundingBox()
+{
+ temp_bb.reset();
+}
+
+void AddToTemporaryBoundingBox(double x, double y, double z)
+{
+ temp_bb += SPoint3(x, y, z);
+ CTX.min[0] = temp_bb.min().x(); CTX.max[0] = temp_bb.max().x();
+ CTX.min[1] = temp_bb.min().y(); CTX.max[1] = temp_bb.max().y();
+ CTX.min[2] = temp_bb.min().z(); CTX.max[2] = temp_bb.max().z();
+ FinishUpBoundingBox();
+}
+
+int ParseFile(char *f, int close, int warn_if_missing)
+{
+ char yyname_old[256], tmp[256];
+ FILE *yyin_old, *fp;
+ int yylineno_old, yyerrorstate_old, numviews_old;
+
+ // add 'b' for pure Windows programs: opening in text mode messes up
+ // fsetpos/fgetpos (used e.g. for user-defined functions)
+ if(!(fp = fopen(f, "rb"))){
+ if(warn_if_missing) Msg(WARNING, "Unable to open file '%s'", f);
+ return 0;
+ }
+
+ strncpy(yyname_old, yyname, 255);
+ yyin_old = yyin;
+ yyerrorstate_old = yyerrorstate;
+ yylineno_old = yylineno;
+ numviews_old = List_Nbr(CTX.post.list);
+
+ strncpy(yyname, f, 255);
+ yyin = fp;
+ yyerrorstate = 0;
+ yylineno = 1;
+
+ fpos_t position;
+ fgetpos(yyin, &position);
+ fgets(tmp, sizeof(tmp), yyin);
+ fsetpos(yyin, &position);
+
+ while(!feof(yyin)){
+ yyparse();
+ if(yyerrorstate > 20){
+ Msg(GERROR, "Too many errors: aborting...");
+ force_yyflush();
+ break;
+ }
+ }
+
+ if(close)
+ fclose(yyin);
+
+ strncpy(yyname, yyname_old, 255);
+ yyin = yyin_old;
+ yyerrorstate = yyerrorstate_old;
+ yylineno = yylineno_old;
+
+ if(List_Nbr(CTX.post.list) != numviews_old){
+#if defined(HAVE_FLTK)
+ UpdateViewsInGUI();
+#endif
+ }
+
+ return 1;
+}
+
+void ParseString(char *str)
+{
+ if(!str) return;
+ FILE *fp;
+ if((fp = fopen(CTX.tmp_filename_fullpath, "w"))) {
+ fprintf(fp, str);
+ fprintf(fp, "\n");
+ fclose(fp);
+ ParseFile(CTX.tmp_filename_fullpath, 1);
+ if(GMODEL) GMODEL->importTHEM();
+ }
+}
+
+void SetProjectName(char *name)
+{
+ char base[356], ext[256];
+ SplitFileName(name, base, ext);
+
+ strncpy(CTX.filename, name, 255);
+ strncpy(CTX.base_filename, base, 255);
+
+#if defined(HAVE_FLTK)
+ if(!CTX.batch) WID->set_title(CTX.filename);
+#endif
+}
+
+int MergeFile(char *name, int warn_if_missing)
+{
+ // FIXME: to be removed
+ // #if defined(HAVE_FOURIER_MODEL)
+ // if(!strcmp(name, "falcon") || !strcmp(name, "ship")){
+ // GMODEL->readFourier(name);
+ // SetBoundingBox();
+ // CTX.mesh.changed = ENT_ALL;
+ // return 1;
+ // }
+ // #endif
+
+ // added 'b' for pure Windows programs, since some of these files
+ // contain binary data
+ FILE *fp = fopen(name, "rb");
+ if(!fp){
+ if(warn_if_missing) Msg(WARNING, "Unable to open file '%s'", name);
+ return 0;
+ }
+
+ char header[256];
+ fgets(header, sizeof(header), fp);
+ fclose(fp);
+
+ Msg(STATUS2, "Reading '%s'", name);
+
+ char ext[256], base[256];
+ SplitFileName(name, base, ext);
+
+#if defined(HAVE_FLTK)
+ if(!CTX.batch) {
+ if(!strcmp(ext, ".gz")) {
+ // the real solution would be to rewrite all our I/O functions in
+ // terms of gzFile, but until then, this is better than nothing
+ if(fl_choice("File '%s' is in gzip format.\n\nDo you want to uncompress it?",
+ "Cancel", "Uncompress", NULL, name)){
+ char tmp[256];
+ sprintf(tmp, "gunzip -c %s > %s", name, base);
+ SystemCall(tmp);
+ if(!strcmp(CTX.filename, name)) // this is the project file
+ SetProjectName(base);
+ return MergeFile(base);
+ }
+ }
+ }
+#endif
+
+ CTX.geom.draw = 0; // don't try to draw the model while reading
+ int status = 0;
+ if(!strcmp(ext, ".stl") || !strcmp(ext, ".STL")){
+ status = GMODEL->readSTL(name, CTX.geom.tolerance);
+ }
+ else if(!strcmp(ext, ".brep") || !strcmp(ext, ".rle") ||
+ !strcmp(ext, ".brp") || !strcmp(ext, ".BRP")){
+ GMODEL->readOCCBREP(std::string(name));
+ }
+ else if(!strcmp(ext, ".iges") || !strcmp(ext, ".IGES") ||
+ !strcmp(ext, ".igs") || !strcmp(ext, ".IGS")){
+ GMODEL->readOCCIGES(std::string(name));
+ }
+ else if(!strcmp(ext, ".step") || !strcmp(ext, ".STEP") ||
+ !strcmp(ext, ".stp") || !strcmp(ext, ".STP")){
+ GMODEL->readOCCSTEP(std::string(name));
+ }
+ else if(!strcmp(ext, ".unv") || !strcmp(ext, ".UNV")){
+ status = GMODEL->readUNV(name);
+ }
+ else if(!strcmp(ext, ".wrl") || !strcmp(ext, ".WRL") ||
+ !strcmp(ext, ".vrml") || !strcmp(ext, ".VRML") ||
+ !strcmp(ext, ".iv") || !strcmp(ext, ".IV")){
+ status = GMODEL->readVRML(name);
+ }
+ else if(!strcmp(ext, ".mesh") || !strcmp(ext, ".MESH")){
+ status = GMODEL->readMESH(name);
+ }
+ else if(!strcmp(ext, ".bdf") || !strcmp(ext, ".BDF") ||
+ !strcmp(ext, ".nas") || !strcmp(ext, ".NAS")){
+ status = GMODEL->readBDF(name);
+ }
+ else if(!strcmp(ext, ".p3d") || !strcmp(ext, ".P3D")){
+ status = GMODEL->readP3D(name);
+ }
+ else if(!strcmp(ext, ".fm") || !strcmp(ext, ".FM")) {
+ status = GMODEL->readF(name);
+ }
+#if defined(HAVE_FLTK)
+ else if(!strcmp(ext, ".pnm") || !strcmp(ext, ".PNM") ||
+ !strcmp(ext, ".pbm") || !strcmp(ext, ".PBM") ||
+ !strcmp(ext, ".pgm") || !strcmp(ext, ".PGM") ||
+ !strcmp(ext, ".ppm") || !strcmp(ext, ".PPM")) {
+ status = read_pnm(name);
+ }
+ else if(!strcmp(ext, ".bmp") || !strcmp(ext, ".BMP")) {
+ status = read_bmp(name);
+ }
+#if defined(HAVE_LIBJPEG)
+ else if(!strcmp(ext, ".jpg") || !strcmp(ext, ".JPG") ||
+ !strcmp(ext, ".jpeg") || !strcmp(ext, ".JPEG")) {
+ status = read_jpeg(name);
+ }
+#endif
+#if defined(HAVE_LIBPNG)
+ else if(!strcmp(ext, ".png") || !strcmp(ext, ".PNG")) {
+ status = read_png(name);
+ }
+#endif
+#endif
+ else {
+ CTX.geom.draw = 1;
+ if(!strncmp(header, "$PTS", 4) || !strncmp(header, "$NO", 3) ||
+ !strncmp(header, "$PARA", 5) || !strncmp(header, "$ELM", 4) ||
+ !strncmp(header, "$MeshFormat", 11)) {
+ status = GMODEL->readMSH(name);
+ }
+ else if(!strncmp(header, "$PostFormat", 11) ||
+ !strncmp(header, "$View", 5)) {
+ status = ReadView(name);
+ }
+ else {
+ status = GMODEL->readGEO(name);
+ }
+ }
+
+ SetBoundingBox();
+
+ CTX.geom.draw = 1;
+ CTX.mesh.changed = ENT_ALL;
+
+ checkHighOrderTriangles ( GMODEL );
+
+ Msg(STATUS2, "Read '%s'", name);
+ return status;
+}
+
+void OpenProject(char *name)
+{
+ if(CTX.threads_lock) {
+ Msg(INFO, "I'm busy! Ask me that later...");
+ return;
+ }
+ CTX.threads_lock = 1;
+
+ GMODEL->destroy();
+ THEM->destroy();
+
+ // Initialize pseudo random mesh generator to the same seed
+ srand(1);
+
+ // temporary hack until we fill GMODEL on the fly during parsing
+ ResetTemporaryBoundingBox();
+
+ SetProjectName(name);
+ MergeFile(name);
+
+ CTX.threads_lock = 0;
+
+#if defined(HAVE_FLTK)
+ if(!CTX.batch) WID->reset_visibility();
+ ZeroHighlight();
+#endif
+}
+
+void OpenProjectMacFinder(const char *filename)
+{
+ static int first = 1;
+ if(first || CTX.batch){
+ // just copy the filename: it will be opened when the GUI is ready
+ // in main()
+ strncpy(CTX.filename, filename, 255);
+ first = 0;
+ }
+ else{
+ OpenProject((char*)filename);
+#if defined(HAVE_FLTK)
+ Draw();
+#endif
+ }
+}
+
diff --git a/configure b/configure
index 22b61a37ca..fe6706e2e5 100755
--- a/configure
+++ b/configure
@@ -866,6 +866,7 @@ Optional Features:
--enable-cgns enable CGNS output (default=no)
--enable-occ enable OpenCascade support (default=no)
--enable-med enable MED support (default=yes)
+ --enable-fm enable support for Fourier models (default=yes)
Optional Packages:
--with-PACKAGE[=ARG] use PACKAGE [ARG=yes]
@@ -1485,6 +1486,11 @@ fi;
if test "${enable_med+set}" = set; then
enableval="$enable_med"
+fi;
+# Check whether --enable-fm or --disable-fm was given.
+if test "${enable_fm+set}" = set; then
+ enableval="$enable_fm"
+
fi;
UNAME=`uname`
@@ -4246,6 +4252,7 @@ fi
fi
fi
+ if test "x$enable_matheval" != "xno"; then
echo "$as_me:$LINENO: checking for ./contrib/MathEval/matheval.cpp" >&5
echo $ECHO_N "checking for ./contrib/MathEval/matheval.cpp... $ECHO_C" >&6
if test "${ac_cv_file___contrib_MathEval_matheval_cpp+set}" = set; then
@@ -4269,43 +4276,43 @@ else
MATHEVAL="no"
fi
- if test "x${MATHEVAL}" = "xyes"; then
- if test "x$enable_matheval" != "xno"; then
- GMSH_DIRS="${GMSH_DIRS} contrib/MathEval"
- GMSH_LIBS="${GMSH_LIBS} -lGmshMathEval"
- FLAGS="-DHAVE_MATH_EVAL ${FLAGS}"
+ if test "x${MATHEVAL}" = "xyes"; then
+ GMSH_DIRS="${GMSH_DIRS} contrib/MathEval"
+ GMSH_LIBS="${GMSH_LIBS} -lGmshMathEval"
+ FLAGS="-DHAVE_MATH_EVAL ${FLAGS}"
fi
fi
- echo "$as_me:$LINENO: checking for ./contrib/FourierModel/model.cpp" >&5
-echo $ECHO_N "checking for ./contrib/FourierModel/model.cpp... $ECHO_C" >&6
-if test "${ac_cv_file___contrib_FourierModel_model_cpp+set}" = set; then
+ if test "x$enable_fm" != "xno"; then
+ echo "$as_me:$LINENO: checking for ./contrib/FourierModel/FM_Face.cpp" >&5
+echo $ECHO_N "checking for ./contrib/FourierModel/FM_Face.cpp... $ECHO_C" >&6
+if test "${ac_cv_file___contrib_FourierModel_FM_Face_cpp+set}" = set; then
echo $ECHO_N "(cached) $ECHO_C" >&6
else
test "$cross_compiling" = yes &&
{ { echo "$as_me:$LINENO: error: cannot check for file existence when cross compiling" >&5
echo "$as_me: error: cannot check for file existence when cross compiling" >&2;}
{ (exit 1); exit 1; }; }
-if test -r "./contrib/FourierModel/model.cpp"; then
- ac_cv_file___contrib_FourierModel_model_cpp=yes
+if test -r "./contrib/FourierModel/FM_Face.cpp"; then
+ ac_cv_file___contrib_FourierModel_FM_Face_cpp=yes
else
- ac_cv_file___contrib_FourierModel_model_cpp=no
+ ac_cv_file___contrib_FourierModel_FM_Face_cpp=no
fi
fi
-echo "$as_me:$LINENO: result: $ac_cv_file___contrib_FourierModel_model_cpp" >&5
-echo "${ECHO_T}$ac_cv_file___contrib_FourierModel_model_cpp" >&6
-if test $ac_cv_file___contrib_FourierModel_model_cpp = yes; then
+echo "$as_me:$LINENO: result: $ac_cv_file___contrib_FourierModel_FM_Face_cpp" >&5
+echo "${ECHO_T}$ac_cv_file___contrib_FourierModel_FM_Face_cpp" >&6
+if test $ac_cv_file___contrib_FourierModel_FM_Face_cpp = yes; then
FOURIER="yes"
else
FOURIER="no"
fi
- if test "x${FOURIER}" = "xyes"; then
- GMSH_DIRS="${GMSH_DIRS} contrib/FourierModel"
- GMSH_LIBS="${GMSH_LIBS} -lGmshFourierModel"
- FLAGS="-DHAVE_FOURIER_MODEL ${FLAGS}"
+ if test "x${FOURIER}" = "xyes"; then
+ GMSH_DIRS="${GMSH_DIRS} contrib/FourierModel"
+ GMSH_LIBS="${GMSH_LIBS} -lGmshFourierModel"
+ FLAGS="-DHAVE_FOURIER_MODEL ${FLAGS}"
+ fi
fi
-
fi
if test "x$enable_gsl" != "xno"; then
diff --git a/configure.in b/configure.in
index c6f6c192c6..6d1190baa4 100644
--- a/configure.in
+++ b/configure.in
@@ -1,4 +1,4 @@
-dnl $Id: configure.in,v 1.124 2007-03-29 16:34:21 geuzaine Exp $
+dnl $Id: configure.in,v 1.125 2007-05-03 08:50:01 geuzaine Exp $
dnl
dnl Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
dnl
@@ -125,6 +125,9 @@ AC_ARG_ENABLE(occ,
AC_ARG_ENABLE(med,
AC_HELP_STRING([--enable-med],
[enable MED support (default=yes)]))
+AC_ARG_ENABLE(fm,
+ AC_HELP_STRING([--enable-fm],
+ [enable support for Fourier models (default=yes)]))
dnl Get the operating system name
UNAME=`uname`
@@ -475,23 +478,24 @@ if test "x$enable_contrib" != "xno"; then
fi
dnl Check for MathEval
- AC_CHECK_FILE(./contrib/MathEval/matheval.cpp, MATHEVAL="yes", MATHEVAL="no")
- if test "x${MATHEVAL}" = "xyes"; then
- if test "x$enable_matheval" != "xno"; then
- GMSH_DIRS="${GMSH_DIRS} contrib/MathEval"
- GMSH_LIBS="${GMSH_LIBS} -lGmshMathEval"
- FLAGS="-DHAVE_MATH_EVAL ${FLAGS}"
+ if test "x$enable_matheval" != "xno"; then
+ AC_CHECK_FILE(./contrib/MathEval/matheval.cpp, MATHEVAL="yes", MATHEVAL="no")
+ if test "x${MATHEVAL}" = "xyes"; then
+ GMSH_DIRS="${GMSH_DIRS} contrib/MathEval"
+ GMSH_LIBS="${GMSH_LIBS} -lGmshMathEval"
+ FLAGS="-DHAVE_MATH_EVAL ${FLAGS}"
fi
fi
dnl Check for FourierModel
- AC_CHECK_FILE(./contrib/FourierModel/model.cpp, FOURIER="yes", FOURIER="no")
- if test "x${FOURIER}" = "xyes"; then
- GMSH_DIRS="${GMSH_DIRS} contrib/FourierModel"
- GMSH_LIBS="${GMSH_LIBS} -lGmshFourierModel"
- FLAGS="-DHAVE_FOURIER_MODEL ${FLAGS}"
+ if test "x$enable_fm" != "xno"; then
+ AC_CHECK_FILE(./contrib/FourierModel/FM_Face.cpp, FOURIER="yes", FOURIER="no")
+ if test "x${FOURIER}" = "xyes"; then
+ GMSH_DIRS="${GMSH_DIRS} contrib/FourierModel"
+ GMSH_LIBS="${GMSH_LIBS} -lGmshFourierModel"
+ FLAGS="-DHAVE_FOURIER_MODEL ${FLAGS}"
+ fi
fi
-
fi
dnl Check for GSL
diff --git a/contrib/FourierModel/ContinuationPatch.cpp b/contrib/FourierModel/ContinuationPatch.cpp
new file mode 100644
index 0000000000..8018987703
--- /dev/null
+++ b/contrib/FourierModel/ContinuationPatch.cpp
@@ -0,0 +1,500 @@
+#include <math.h>
+#include "Message.h"
+#include "ContinuationPatch.h"
+
+ContinuationPatch::ContinuationPatch(PatchInfo *PI, int derivative) : Patch(PI)
+{
+ _NcrossT[0] = _PI->normal[1] * _PI->tangent[2] -
+ _PI->normal[2] * _PI->tangent[1];
+ _NcrossT[1] = _PI->normal[2] * _PI->tangent[0] -
+ _PI->normal[0] * _PI->tangent[2];
+ _NcrossT[2] = _PI->normal[0] * _PI->tangent[1] -
+ _PI->normal[1] * _PI->tangent[0];
+
+ _derivative = derivative;
+
+ _uModes = _PI->nModes[0];
+ _vModes = _PI->nModes[1];
+
+ _periodU = (1-_PI->periodic[0]) + 1;
+ _periodV = (1-_PI->periodic[1]) + 1;
+
+ //Msg::Info("Period : %g %g",_periodU,_periodV);
+
+ if ((_uModes % 2) == 0) {
+ _uModesLower = - _uModes/2;
+ _uModesUpper = _uModes/2 - 1;
+ }
+ else {
+ _uModesLower = - (_uModes - 1)/2;
+ _uModesUpper = (_uModes - 1)/2;
+ }
+ if ((_vModes % 2) == 0) {
+ _vModesLower = - _vModes/2;
+ _vModesUpper = _vModes/2 - 1;
+ }
+ else {
+ _vModesLower = - (_vModes - 1)/2;
+ _vModesUpper = (_vModes - 1)/2;
+ }
+
+ // Initialize Data
+ _coeffData = new std::complex<double>*[_uModes];
+ for(int j = 0; j < _uModes; j++){
+ _coeffData[j] = new std::complex<double>[_vModes];
+ for(int k = 0; k < _vModes; k++)
+ _coeffData[j][k] = _PI->coeff[j][k];
+ }
+
+ // Check if we need to interpolate the derivative(s)
+ if(_derivative){
+ // Initialize _fineDeriv and _fineDeriv2 to zero
+ if(_derivative & 1){
+ _coeffDerivU = new std::complex<double>*[_uModes];
+ _coeffDerivV = new std::complex<double>*[_uModes];
+ for(int j = 0; j < _uModes; j++){
+ _coeffDerivU[j] = new std::complex<double>[_vModes];
+ _coeffDerivV[j] = new std::complex<double>[_vModes];
+ for(int k = 0; k < _vModes; k++){
+ _coeffDerivU[j][k] = 0.;
+ _coeffDerivV[j][k] = 0.;
+ }
+ }
+ }
+
+ if(_derivative & 2){
+ _coeffDerivUU = new std::complex<double>*[_uModes];
+ _coeffDerivVV = new std::complex<double>*[_uModes];
+ _coeffDerivUV = new std::complex<double>*[_uModes];
+ for(int j = 0; j < _uModes; j++){
+ _coeffDerivUU[j] = new std::complex<double>[_vModes];
+ _coeffDerivVV[j] = new std::complex<double>[_vModes];
+ _coeffDerivUV[j] = new std::complex<double>[_vModes];
+ for(int k = 0; k < _vModes; k++){
+ _coeffDerivUU[j][k] = 0.;
+ _coeffDerivVV[j][k] = 0.;
+ _coeffDerivUV[j][k] = 0.;
+ }
+ }
+ }
+
+ // Copy the Fourier coefficients into _coeffDeriv and _coeffDeriv2
+ std::complex<double> I(0., 1.);
+ for(int j = 0; j < _uModes; j++){
+ for(int k = 0; k < _vModes; k++){
+ int J = j+_uModesLower;
+ int K = k+_vModesLower;
+ if(_derivative & 1){
+ _coeffDerivU[j][k] = (2 * M_PI * J * I / _periodU) *_coeffData[j][k];
+ _coeffDerivV[j][k] = (2 * M_PI * K * I / _periodV) *_coeffData[j][k];
+ }
+ if(_derivative & 2){
+ _coeffDerivUU[j][k] = - (4 * M_PI * M_PI * J * J /
+ (_periodU * _periodU)) * _coeffData[j][k];
+ _coeffDerivVV[j][k] = - (4 * M_PI * M_PI * K * K /
+ (_periodV * _periodV)) * _coeffData[j][k];
+ _coeffDerivUV[j][k] = - (4 * M_PI * M_PI * J * K /
+ (_periodU * _periodV)) * _coeffData[j][k];
+ }
+ }
+ }
+ }
+ // Initialize interpolation variables
+ _tmpCoeff = std::vector< std::complex<double> >(_vModes);
+ _tmpInterp = std::vector< std::complex<double> >(_uModes);
+}
+
+ContinuationPatch::~ContinuationPatch()
+{
+ for(int j = 0; j < _uModes; j++)
+ delete [] _coeffData[j];
+ delete [] _coeffData;
+
+ if(_coeffDerivU){
+ for(int j = 0; j < _uModes; j++)
+ delete [] _coeffDerivU[j];
+ delete [] _coeffDerivU;
+ }
+ if(_coeffDerivV){
+ for(int j = 0; j < _uModes; j++)
+ delete [] _coeffDerivV[j];
+ delete [] _coeffDerivV;
+ }
+ if(_coeffDerivUU){
+ for(int j = 0; j < _uModes; j++)
+ delete [] _coeffDerivUU[j];
+ delete [] _coeffDerivUU;
+ }
+ if(_coeffDerivVV){
+ for(int j = 0; j < _uModes; j++)
+ delete [] _coeffDerivVV[j];
+ delete [] _coeffDerivVV;
+ }
+ if(_coeffDerivUV){
+ for(int j = 0; j < _uModes; j++)
+ delete [] _coeffDerivUV[j];
+ delete [] _coeffDerivUV;
+ }
+}
+
+std::complex<double> ContinuationPatch::
+_PolyEval(std::vector< std::complex<double> > _coeff, std::complex<double> x)
+{
+ int _polyOrder = _coeff.size()-1;
+ std::complex<double> out = 0.;
+
+ out = x * _coeff[_polyOrder];
+ for (int i = _polyOrder - 1; i > 0; i--)
+ out = x * (out + _coeff[i]);
+ out = out + _coeff[0];
+
+ return out;
+}
+
+std::complex<double> ContinuationPatch::
+ _Interpolate(double u, double v, int uDer, int vDer)
+{
+ //Msg::Info("%d %d %d",uDer,vDer,_derivative);
+ if (((uDer==2 || vDer==2 || (uDer==1 && vDer==1)) && !(_derivative & 2) ) ||
+ ((uDer==1 || vDer==1) && !(_derivative & 1)) ||
+ (uDer<0 || uDer>2 || vDer<0 || vDer>2) ) {
+ Msg::Error("Derivative data not available: check contructor call %d %d %d",
+ uDer,vDer,_derivative);
+ return 0.;
+ }
+
+ double epsilon = 1e-12;
+ if(u < 0. - epsilon || u > 1. + epsilon){
+ Msg::Error("Trying to interpolate outside interval: (u,v)=(%.16g,%.16g) "
+ "not in [%g,%g]x[%g,%g]", u, v, 0., 1., 0., 1.);
+ }
+
+ // Interpolate to find value at (u,v)
+ for(int j = 0; j < _uModes; j++){
+ for(int k = 0; k < _vModes; k++){
+ std::complex<double> tmp;
+ if(uDer == 0 && vDer == 0)
+ tmp = _coeffData[j][k];
+ else if(uDer == 1 && vDer == 0)
+ tmp = _coeffDerivU[j][k];
+ else if(uDer == 0 && vDer == 1)
+ tmp = _coeffDerivV[j][k];
+ else if(uDer == 2 && vDer == 0)
+ tmp = _coeffDerivUU[j][k];
+ else if(uDer == 0 && vDer == 2)
+ tmp = _coeffDerivVV[j][k];
+ else
+ tmp = _coeffDerivUV[j][k];
+ _tmpCoeff[k] = tmp;
+ }
+ std::complex<double> y(cos(2 * M_PI * v / _periodV),
+ sin(2 * M_PI * v / _periodV));
+ _tmpInterp[j] = _PolyEval(_tmpCoeff, y);
+ _tmpInterp[j] *= std::complex<double>
+ (cos(2 * M_PI * _vModesLower * v / _periodV),
+ sin(2 * M_PI * _vModesLower * v / _periodV));
+ }
+ std::complex<double> x(cos(2 * M_PI * u / _periodU),
+ sin(2 * M_PI * u / _periodU));
+ return _PolyEval(_tmpInterp, x) * std::complex<double>
+ (cos(2 * M_PI * _uModesLower * u / _periodU),
+ sin(2 * M_PI * _uModesLower * u / _periodU));
+}
+
+void ContinuationPatch::P(double u, double v, double &x, double &y, double &z,
+ double &nx, double &ny, double &nz)
+{
+ if (!strcmp(_PI->projection,"cylinder")) {
+ x = _PI->origin[0] + _PI->normal[0] * _PI->scale[1] * v;
+ y = _PI->origin[1] + _PI->normal[1] * _PI->scale[1] * v;
+ z = _PI->origin[2] + _PI->normal[2] * _PI->scale[1] * v;
+
+ nx = _PI->tangent[0] * cos(_PI->scale[0] * (u - 0.5)) +
+ _NcrossT[0] * sin(_PI->scale[0] * (u - 0.5));
+ ny = _PI->tangent[1] * cos(_PI->scale[0] * (u - 0.5)) +
+ _NcrossT[1] * sin(_PI->scale[0] * (u - 0.5));
+ nz = _PI->tangent[2] * cos(_PI->scale[0] * (u - 0.5)) +
+ _NcrossT[2] * sin(_PI->scale[0] * (u - 0.5));
+
+ double norm = sqrt(nx * nx + ny * ny + nz * nz);
+ nx /= norm;
+ ny /= norm;
+ nz /= norm;
+ }
+ else if (!strcmp(_PI->projection,"cylinder")) {
+ x = y = z = nx = ny = nz = 0;
+ }
+ else {
+ Msg::Info("Unknown projection surface");
+ x = y = z = nx = ny = nz = 0;
+ }
+}
+
+bool ContinuationPatch::pInverse(double x,double y,double z,
+ double &u,double &v)
+{
+ if (!strcmp(_PI->projection,"cylinder")) {
+ double t = (x - _PI->origin[0]) * _PI->normal[0] +
+ (y - _PI->origin[1]) * _PI->normal[1] +
+ (z - _PI->origin[2]) * _PI->normal[2];
+ v = t / _PI->scale[1];
+ double n[3];
+ n[0] = x - (_PI->origin[0] + t * _PI->normal[0]);
+ n[1] = y - (_PI->origin[1] + t * _PI->normal[1]);
+ n[2] = z - (_PI->origin[2] + t * _PI->normal[2]);
+ double norm = sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
+ for (int i=0;i<3;i++)
+ n[i] /= norm;
+ u = atan2((n[0]*_NcrossT[0]+n[1]*_NcrossT[1]+n[2]*_NcrossT[2]),
+ (n[0]*_PI->tangent[0]+n[1]*_PI->tangent[1]+
+ n[2]*_PI->tangent[2]));
+ u /= _PI->scale[0];
+ u += 0.5;
+
+ }
+ else if (!strcmp(_PI->projection,"cylinder")) {
+ u = v = 0.;
+ }
+ else {
+ Msg::Info("Unknown projection surface");
+ u = v = 0.;
+ }
+}
+
+void ContinuationPatch::Dpdu(double u, double v,
+ double &x, double &y, double &z,
+ double &nx, double &ny, double &nz)
+{
+ if (!strcmp(_PI->projection,"cylinder")) {
+ x = y = z = 0.;
+
+ double n[3], nu[3];
+ for (int i=0;i<3;i++) {
+ n[i] = _PI->tangent[i] * cos(_PI->scale[0] * u) +
+ _NcrossT[i] * sin(_PI->scale[0] * u);
+ nu[i] = _PI->scale[0] * (- _PI->tangent[i] * sin(_PI->scale[0] * u) +
+ _NcrossT[i] * cos(_PI->scale[0] * u));
+ }
+ double norm = sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
+ double NdotNu = n[0] * nu[0] + n[1] * nu[1] + n[2] * nu[2];
+ nx = nu[0] / norm - NdotNu * n[0] / (norm * norm * norm);
+ ny = nu[1] / norm - NdotNu * n[1] / (norm * norm * norm);
+ nz = nu[2] / norm - NdotNu * n[2] / (norm * norm * norm);
+ }
+ else if (!strcmp(_PI->projection,"cylinder")) {
+ x = y = z = nx = ny = nz = 0;
+ }
+ else {
+ Msg::Info("Unknown projection surface");
+ x = y = z = nx = ny = nz = 0;
+ }
+}
+
+void ContinuationPatch::Dpdv(double u, double v,
+ double &x, double &y, double &z,
+ double &nx, double &ny, double &nz)
+{
+ if (!strcmp(_PI->projection,"cylinder")) {
+ x = _PI->normal[0] * _PI->scale[1];
+ y = _PI->normal[1] * _PI->scale[1];
+ z = _PI->normal[2] * _PI->scale[1];
+
+ nx = ny = nz = 0.;
+ }
+ else if (!strcmp(_PI->projection,"cylinder")) {
+ x = y = z = nx = ny = nz = 0;
+ }
+ else {
+ Msg::Info("Unknown projection surface");
+ x = y = z = nx = ny = nz = 0;
+ }
+}
+
+void ContinuationPatch::Dpdpdudu(double u,double v,
+ double &x, double &y, double &z,
+ double &nx, double &ny, double &nz)
+{
+ if (!strcmp(_PI->projection,"cylinder")) {
+ x = y = z = 0.;
+ double n[3], nu[3], nuu[3];
+ for (int i=0;i<3;i++) {
+ n[i] = _PI->tangent[i] * cos(_PI->scale[0] * (u-0.5)) +
+ _NcrossT[i] * sin(_PI->scale[0] * (u-0.5));
+ nu[i] = _PI->scale[0] *
+ (- _PI->tangent[i] * sin(_PI->scale[0] * (u-0.5)) +
+ _NcrossT[i] * cos(_PI->scale[0] * (u-0.5)));
+ nuu[i] = - _PI->scale[0] * _PI->scale[0] *
+ (_PI->tangent[i] * cos(_PI->scale[0] * (u-0.5)) +
+ _NcrossT[i] * sin(_PI->scale[0] * (u-0.5)));
+ }
+ double norm = sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
+ double NdotNu = n[0] * nu[0] + n[1] * nu[1] + n[2] * nu[2];
+ double NudotNu = nu[0] * nu[0] + nu[1] * nu[1] + nu[2] * nu[2];
+ double NdotNuu = n[0] * nuu[0] + n[1] * nuu[1] + n[2] * nuu[2];
+
+ nx = nuu[0] / norm - (NdotNu * nu[0] + (NdotNuu + NudotNu) * n[0] +
+ NdotNu * nu[0]) / (norm * norm * norm) +
+ 3 * NdotNu * NdotNu * n[0] / (norm * norm * norm * norm * norm);
+ ny = nuu[1] / norm - (NdotNu * nu[1] + (NdotNuu + NudotNu) * n[1] +
+ NdotNu * nu[1]) / (norm * norm * norm) +
+ 3 * NdotNu * NdotNu * n[1] / (norm * norm * norm * norm * norm);
+ nz = nuu[2] / norm - (NdotNu * nu[2] + (NdotNuu + NudotNu) * n[2] +
+ NdotNu * nu[2]) / (norm * norm * norm) +
+ 3 * NdotNu * NdotNu * n[2] / (norm * norm * norm * norm * norm);
+ }
+ else if (!strcmp(_PI->projection,"cylinder")) {
+ x = y = z = nx = ny = nz = 0;
+ }
+ else {
+ Msg::Info("Unknown projection surface");
+ x = y = z = nx = ny = nz = 0;
+ }
+}
+
+void ContinuationPatch::Dpdpdudv(double u, double v,
+ double &x, double &y, double &z,
+ double &nx, double &ny, double &nz)
+{
+ if (!strcmp(_PI->projection,"cylinder")) {
+ x = y = z = nx = ny = nz = 0.;
+ }
+ else if (!strcmp(_PI->projection,"cylinder")) {
+ x = y = z = nx = ny = nz = 0;
+ }
+ else {
+ Msg::Info("Unknown projection surface");
+ x = y = z = nx = ny = nz = 0;
+ }
+}
+
+void ContinuationPatch::Dpdpdvdv(double u, double v,
+ double &x, double &y, double &z,
+ double &nx, double &ny, double &nz)
+{
+ if (!strcmp(_PI->projection,"cylinder")) {
+ x = y = z = nx = ny = nz = 0.;
+ }
+ else if (!strcmp(_PI->projection,"cylinder")) {
+ x = y = z = nx = ny = nz = 0;
+ }
+ else {
+ Msg::Info("Unknown projection surface");
+ x = y = z = nx = ny = nz = 0;
+ }
+}
+
+void ContinuationPatch::
+F(double u, double v, double &x, double &y, double &z)
+{
+ double px, py, pz, nx, ny, nz, d;
+ P(u,v,px,py,pz,nx,ny,nz);
+ d = _Interpolate(u, v, 0, 0).real();
+ x = px + d * nx;
+ y = py + d * ny;
+ z = pz + d * nz;
+}
+
+bool ContinuationPatch::
+Inverse(double x,double y,double z,double &u,double &v)
+{
+ pInverse(x,y,z,u,v);
+}
+
+void ContinuationPatch::
+Dfdu(double u, double v, double &x, double &y, double &z)
+{
+ double px, py, pz, nx, ny, nz, d;
+ double pxu, pyu, pzu, nxu, nyu, nzu, du;
+ P(u,v,px,py,pz,nx,ny,nz);
+ Dpdu(u,v,pxu,pyu,pzu,nxu,nyu,nzu);
+ d = _Interpolate(u, v, 0, 0).real();
+ du = _Interpolate(u, v, 1, 0).real();
+
+ x = pxu + du * nx + d * nxu;
+ y = pyu + du * ny + d * nyu;
+ z = pzu + du * nz + d * nzu;
+}
+
+void ContinuationPatch::
+Dfdv(double u, double v, double &x, double &y, double &z)
+{
+ double px, py, pz, nx, ny, nz, d;
+ double pxv, pyv, pzv, nxv, nyv, nzv, dv;
+ P(u,v,px,py,pz,nx,ny,nz);
+ Dpdu(u,v,pxv,pyv,pzv,nxv,nyv,nzv);
+ d = _Interpolate(u, v, 0, 0).real();
+ dv = _Interpolate(u, v, 0, 1).real();
+
+ x = pxv + dv * nx + d * nxv;
+ y = pyv + dv * ny + d * nyv;
+ z = pzv + dv * nz + d * nzv;
+}
+
+void ContinuationPatch::
+Dfdfdudu(double u, double v, double &x, double &y, double &z)
+{
+ double px, py, pz, nx, ny, nz, d;
+ double pxu, pyu, pzu, nxu, nyu, nzu, du;
+ double pxuu, pyuu, pzuu, nxuu, nyuu, nzuu, duu;
+
+ P(u,v,px,py,pz,nx,ny,nz);
+ Dpdu(u,v,pxu,pyu,pzu,nxu,nyu,nzu);
+ Dpdpdudu(u,v,pxuu,pyuu,pzuu,nxuu,nyuu,nzuu);
+
+ d = _Interpolate(u, v, 0, 0).real();
+ du = _Interpolate(u, v, 1, 0).real();
+ duu = _Interpolate(u, v, 2, 0).real();
+
+ x = pxuu + duu * nx + du * nxu + du * nxu + d * nxuu;
+ y = pyuu + duu * ny + du * nyu + du * nyu + d * nyuu;
+ z = pzuu + duu * nz + du * nzu + du * nzu + d * nzuu;
+}
+
+void ContinuationPatch::
+Dfdfdvdv(double u, double v, double &x, double &y, double &z)
+{
+ double px, py, pz, nx, ny, nz, d;
+ double pxv, pyv, pzv, nxv, nyv, nzv, dv;
+ double pxvv, pyvv, pzvv, nxvv, nyvv, nzvv, dvv;
+
+ P(u,v,px,py,pz,nx,ny,nz);
+ Dpdv(u,v,pxv,pyv,pzv,nxv,nyv,nzv);
+ Dpdpdvdv(u,v,pxvv,pyvv,pzvv,nxvv,nyvv,nzvv);
+
+ d = _Interpolate(u, v, 0, 0).real();
+ dv = _Interpolate(u, v, 0, 1).real();
+ dvv = _Interpolate(u, v, 0, 2).real();
+
+ x = pxvv + dvv * nx + dv * nxv + dv * nxv + d * nxvv;
+ y = pyvv + dvv * ny + dv * nyv + dv * nyv + d * nyvv;
+ z = pzvv + dvv * nz + dv * nzv + dv * nzv + d * nzvv;
+}
+
+void ContinuationPatch::
+Dfdfdudv(double u, double v, double &x, double &y, double &z)
+{
+ double px, py, pz, nx, ny, nz, d;
+ double pxu, pyu, pzu, nxu, nyu, nzu, du;
+ double pxv, pyv, pzv, nxv, nyv, nzv, dv;
+ double pxuv, pyuv, pzuv, nxuv, nyuv, nzuv, duv;
+
+ P(u,v,px,py,pz,nx,ny,nz);
+ Dpdu(u,v,pxu,pyu,pzu,nxu,nyu,nzu);
+ Dpdv(u,v,pxv,pyv,pzv,nxv,nyv,nzv);
+ Dpdpdudv(u,v,pxuv,pyuv,pzuv,nxuv,nyuv,nzuv);
+
+ d = _Interpolate(u, v, 0, 0).real();
+ du = _Interpolate(u, v, 1, 0).real();
+ dv = _Interpolate(u, v, 1, 0).real();
+ duv = _Interpolate(u, v, 1, 1).real();
+
+ x = pxuv + duv * nx + du * nxv + dv * nxu + d * nxuv;
+ y = pyuv + duv * ny + du * nyv + dv * nyu + d * nyuv;
+ z = pzuv + duv * nz + du * nzv + dv * nzu + d * nzuv;
+}
+
+double ContinuationPatch::GetPou(double u, double v)
+{
+ return 1.;
+}
+
diff --git a/contrib/FourierModel/ContinuationPatch.h b/contrib/FourierModel/ContinuationPatch.h
new file mode 100644
index 0000000000..4e5c31f84b
--- /dev/null
+++ b/contrib/FourierModel/ContinuationPatch.h
@@ -0,0 +1,63 @@
+#ifndef _CONTUNATION_PATCH_H_
+#define _CONTINUATION_PATCH_H_
+
+#include "Patch.h"
+
+// The base class for the patches
+class ContinuationPatch : public Patch {
+ private:
+ // third axis
+ double _NcrossT[3];
+ // bitfield telling if we also interpolate the derivative(s)
+ int _derivative;
+ // Number of Fourier Modes
+ int _uModes, _vModes;
+ // Period Information
+ double _periodU, _periodV;
+ // Limits of Fourier Series
+ int _uModesLower, _uModesUpper, _vModesLower, _vModesUpper;
+ // data (and its first 2 derivatives)
+ std::complex<double> **_coeffData, **_coeffDerivU, **_coeffDerivV;
+ std::complex<double> **_coeffDerivUU, **_coeffDerivVV, **_coeffDerivUV;
+ // temporary interpolation variables
+ std::vector< std::complex<double> > _tmpCoeff, _tmpInterp;
+ // polynomial evaluator
+ std::complex<double> _PolyEval(std::vector< std::complex<double> > _coeff,
+ std::complex<double> x);
+ // interpolation wrapper
+ std::complex<double> _Interpolate(double u,double v,int uDer=0,int vDer=0);
+ void P(double u, double v, double &x, double &y, double &z,
+ double &nx, double &ny, double &nz);
+ bool pInverse(double x,double y,double z,double &u,double &v);
+ void Dpdu(double u, double v, double &x, double &y, double &z,
+ double &nx, double &ny, double &nz);
+ void Dpdv(double u, double v, double &x, double &y, double &z,
+ double &nx, double &ny, double &nz);
+ void Dpdpdudu(double u,double v,double &x,double &y,double &z,
+ double &nx, double &ny, double &nz);
+ void Dpdpdudv(double u,double v,double &x,double &y,double &z,
+ double &nx, double &ny, double &nz);
+ void Dpdpdvdv(double u,double v,double &x,double &y,double &z,
+ double &nx, double &ny, double &nz);
+ public:
+ ContinuationPatch(PatchInfo* PI, int derivative = 0);
+ virtual ~ContinuationPatch();
+
+ // These are the virtual functions that must be provided by all
+ // derived patches: GetPou() returns the original smooth
+ // (non-normalized) cutoff on the patch; F() and Inverse() implement
+ // the mapping f: (u,v)->(x,y,z) and its inverse; and the Df*() and Dn*()
+ // functions return the derivatives of the mapping f and unit normal n
+ // with respect to u and v
+ virtual double GetPou(double u, double v);
+ virtual void F(double u, double v, double &x, double &y, double &z);
+ virtual bool Inverse(double x,double y,double z,double &u,double &v);
+ virtual void Dfdu(double u, double v, double &x, double &y, double &z);
+ virtual void Dfdv(double u, double v, double &x, double &y, double &z);
+ virtual void Dfdfdudu(double u,double v,double &x,double &y,double &z);
+ virtual void Dfdfdudv(double u,double v,double &x,double &y,double &z);
+ virtual void Dfdfdvdv(double u,double v,double &x,double &y,double &z);
+
+};
+
+#endif
diff --git a/contrib/FourierModel/Curve.cpp b/contrib/FourierModel/Curve.cpp
new file mode 100644
index 0000000000..81191ac27f
--- /dev/null
+++ b/contrib/FourierModel/Curve.cpp
@@ -0,0 +1,242 @@
+#include <math.h>
+#include "Message.h"
+#include "Curve.h"
+
+Curve::Curve(IntersectionInfo* II, std::vector<Patch*> patches) : _II(II)
+{
+ _h = 1.e-6;
+ _tol = 1.e-12;
+
+ for (int i=0;i<patches.size();i++) {
+ if (patches[i]->GetTag() == _II->intersectingPatches[0].patchTag) {
+ _patch0 = patches[i];
+ break;
+ }
+ }
+ for (int i=0;i<patches.size();i++)
+ if (patches[i]->GetTag() == _II->intersectingPatches[1].patchTag) {
+ _patch1 = patches[i];
+ break;
+ }
+ /*
+ int n1=64;
+ int n2=64;
+ double h1 = 1000./(n1-1);
+ double h2 = 1000./(n2-1);
+ for (int j=0;j<n1;j++) {
+ for (int k=0;k<n2;k++) {
+ double u = j*h1;
+ double v = k*h2;
+ double x, y, z;
+ _patch1->F(u,v,x,y,z);
+ printf("x(%d,%d) = %g; y(%d,%d) = %g; z(%d,%d) = %g;\n",
+ j+1,k+1,x,j+1,k+1,y,j+1,k+1,z);
+ }
+ }
+ */
+ double u,v,x,y,z;
+ x = _II->SP[0]; y = _II->SP[1]; z = _II->SP[2];
+ _patch0->Inverse(x,y,z,u,v);
+ _SP0[0] = u;
+ _SP0[1] = v;
+ _patch1->Inverse(x,y,z,u,v);
+ _SP1[0] = u;
+ _SP1[1] = v;
+ //printf("SPx = %g; SPy = %g; SPz = %g;\n",x,y,z);
+
+ x = _II->EP[0]; y = _II->EP[1]; z = _II->EP[2];
+ _patch0->Inverse(x,y,z,u,v);
+ _EP0[0] = u;
+ _EP0[1] = v;
+ _patch1->Inverse(x,y,z,u,v);
+ _EP1[0] = u;
+ _EP1[1] = v;
+ //printf("EPx = %g; EPy = %g; EPz = %g;\n",x,y,z);
+}
+
+void Curve::F(double t, double &x, double &y, double &z)
+{
+ double u0, v0, u1, v1;
+ //Msg::Info("patch0 : %d",_patch0->GetTag());
+ //Msg::Info("patch1 : %d",_patch1->GetTag());
+ //Msg::Info("%g,%g,%g",_II->SP[0],_II->SP[1],_II->SP[2]);
+ //Msg::Info("%g,%g,%g",_II->EP[0],_II->EP[1],_II->EP[2]);
+ //Msg::Info("%g,%g,%g,%g",_SP0[0],_SP0[1],_EP0[0],_EP0[1]);
+ //Msg::Info("%g,%g,%g,%g",_SP1[0],_SP1[1],_EP1[0],_EP1[1]);
+
+ if ((std::abs(_SP0[0]-_EP0[0])<1.e-12) &&
+ (_patch0->_PI->periodic[0] == 1)) {
+ u0 = _SP0[0] + t * (1. + _EP0[0] - _SP0[0]);
+ if (u0 > 1.)
+ u0 -= floor(u0);
+ }
+ else
+ u0 = _SP0[0] + t * (_EP0[0] - _SP0[0]);
+
+ if ((std::abs(_SP0[1]-_EP0[1])<1.e-12) &&
+ (_patch0->_PI->periodic[1] == 1)) {
+ v0 = _SP0[1] + t * (1. + _EP0[1] - _SP0[1]);
+ if (u0 > 1.)
+ v0 -= floor(v0);
+ }
+ else
+ v0 = _SP0[1] + t * (_EP0[1] - _SP0[1]);
+
+ if ((std::abs(_SP1[0]-_EP1[0])<1.e-12) &&
+ (_patch1->_PI->periodic[0] == 1)) {
+ u1 = _SP1[0] + t * (1. + _EP1[0] - _SP1[0]);
+ if (u1 > 1.)
+ u1 -= floor(u1);
+ }
+ else {
+ u1 = _SP1[0] + t * (_EP1[0] - _SP1[0]);
+ }
+
+ if ((std::abs(_SP1[1]-_EP1[1])<1.e-12) &&
+ (_patch1->_PI->periodic[1] == 1)) {
+ v1 = _SP1[1] + t * (1. + _EP1[1] - _SP1[1]);
+ if (u1 > 1.)
+ v1 -= floor(v1);
+ }
+ else
+ v1 = _SP1[1] + t * (_EP1[1] - _SP1[1]);
+
+ //Msg::Info("%g,%g,%g,%g",u0,v0,u1,v1);
+
+ double x0, y0, z0;
+ _patch0->F(u0,v0,x0,y0,z0);
+
+ //Msg::Info("%g,%g,%g",x0,y0,z0);
+
+ double x1, y1, z1;
+ _patch1->F(u1,v1,x1,y1,z1);
+
+ //Msg::Info("%g,%g,%g",x1,y1,z1);
+
+ double r,u,v,rPlus,uPlus,vPlus;
+
+ r = u0;
+ u = u1;
+ v = v1;
+
+ double F[3];
+ F[0] = x1 - x0;
+ F[1] = y1 - y0;
+ F[2] = z1 - z0;
+
+ while ((std::abs(F[0])>_tol) || (std::abs(F[1])>_tol) ||
+ (std::abs(F[2])>_tol)) {
+ rPlus = r + _h;
+ uPlus = u + _h;
+ vPlus = v + _h;
+ if (_patch0->_PI->periodic[0] == 1)
+ rPlus -= std::floor(rPlus);
+ if (_patch1->_PI->periodic[0] == 1)
+ uPlus -= std::floor(uPlus);
+ if (_patch1->_PI->periodic[1] == 1)
+ vPlus -= floor(vPlus);
+ //Msg::Info(" F : %g,%g,%g",F[0],F[1],F[2]);
+ //Msg::Info(" R : %g,%g,%g",x0,y0,z0);
+ double x0rPlus, y0rPlus, z0rPlus;
+ _patch0->F(rPlus,v0,x0rPlus,y0rPlus,z0rPlus);
+ double x1uPlus, y1uPlus, z1uPlus;
+ _patch1->F(uPlus,v,x1uPlus,y1uPlus,z1uPlus);
+ double x1vPlus, y1vPlus, z1vPlus;
+ _patch1->F(u,vPlus,x1vPlus,y1vPlus,z1vPlus);
+
+ double Df[3][3];
+ Df[0][0] = - (x0rPlus - x0) / _h;
+ Df[0][1] = (x1uPlus - x1) / _h;
+ Df[0][2] = (x1vPlus - x1) / _h;
+ Df[1][0] = - (y0rPlus - y0) / _h;
+ Df[1][1] = (y1uPlus - y1) / _h;
+ Df[1][2] = (y1vPlus - y1) / _h;
+ Df[2][0] = - (z0rPlus - z0) / _h;
+ Df[2][1] = (z1uPlus - z1) / _h;
+ Df[2][2] = (z1vPlus - z1) / _h;
+
+ //Msg::Info("D1: %g,%g,%g",Df[0][0],Df[0][1],Df[0][2]);
+ //Msg::Info("D2: %g,%g,%g",Df[1][0],Df[1][1],Df[1][2]);
+ //Msg::Info("D3: %g,%g,%g",Df[2][0],Df[2][1],Df[2][2]);
+
+ double det =
+ Df[0][0] * (Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1]) +
+ Df[0][1] * (Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2]) +
+ Df[0][2] * (Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0]);
+
+ //Msg::Info("det: %g",det);
+
+ double update[3];
+ update[0] =
+ (Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1]) * F[0] +
+ (Df[0][2] * Df[2][1] - Df[0][1] * Df[2][2]) * F[1] +
+ (Df[0][1] * Df[1][2] - Df[0][2] * Df[1][1]) * F[2];
+ update[1] =
+ (Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2]) * F[0] +
+ (Df[0][0] * Df[2][2] - Df[0][2] * Df[2][0]) * F[1] +
+ (Df[0][2] * Df[1][0] - Df[0][0] * Df[1][2]) * F[2];
+ update[2] =
+ (Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0]) * F[0] +
+ (Df[0][1] * Df[2][0] - Df[0][0] * Df[2][1]) * F[1] +
+ (Df[0][0] * Df[1][1] - Df[0][1] * Df[1][0]) * F[2];
+
+ //Msg::Info("U: %g,%g,%g",update[0],update[1],update[2]);
+
+ r -= update[0] / det;
+ u -= update[1] / det;
+ v -= update[2] / det;
+
+ //Msg::Info("U/det: %g,%g,%g",update[0]/det,update[1]/det,update[2]/det);
+
+ //Msg::Info("uv: %g,%g,%g",r,u,v);
+
+ if (_patch0->_PI->periodic[0] == 1)
+ r -= std::floor(r);
+ if (_patch1->_PI->periodic[0] == 1)
+ u -= std::floor(u);
+ if (_patch1->_PI->periodic[1] == 1)
+ v -= floor(v);
+
+ //Msg::Info("UV : %g,%g,%g",r,u,v);
+
+ //rPlus = r + _h;
+ //uPlus = u + _h;
+ //vPlus = v + _h;
+
+ _patch0->F(r,v0,x0,y0,z0);
+ _patch1->F(u,v,x1,y1,z1);
+
+ F[0] = x1 - x0;
+ F[1] = y1 - y0;
+ F[2] = z1 - z0;
+ }
+ x = x0; y = y0; z = z0;
+}
+
+bool Curve::Inverse(double x,double y,double z,double &t)
+{
+ double u0, v0;
+ _patch0->Inverse(x,y,z,u0,v0);
+
+ if ((std::abs(_SP0[1]-_EP0[1])<1.e-12) &&
+ (_patch0->_PI->periodic[1] == 1)) {
+ t = (v0 - _SP0[1] > 0 ? v0 - _SP0[1] : 1. + v0 - _SP0[1]);
+ }
+ else
+ t = (v0 - _SP0[1]) / (_EP0[1] - _SP0[1]);
+
+ if ((t < 0.) || (t > 1.))
+ return false;
+ else
+ return true;
+}
+
+double Curve::GetPou(double t)
+{
+ return 1.;
+}
+
+int Curve::GetTag()
+{
+ return _II->tag;
+}
diff --git a/contrib/FourierModel/Curve.h b/contrib/FourierModel/Curve.h
new file mode 100644
index 0000000000..8606061317
--- /dev/null
+++ b/contrib/FourierModel/Curve.h
@@ -0,0 +1,37 @@
+#ifndef _CURVE_H_
+#define _CURVE_H_
+
+#include "Patch.h"
+#include "FM_Info.h"
+
+// The base class for the patches
+class Curve {
+ private:
+ double _h, _tol;
+ protected:
+ // Patches
+ Patch* _patch0;
+ Patch* _patch1;
+ // End Points
+ double _SP0[2],_EP0[2];
+ double _SP1[2],_EP1[2];
+
+ public:
+ // Intersection Information
+ IntersectionInfo* _II;
+ Curve(IntersectionInfo* II, std::vector<Patch*> patches);
+ virtual ~Curve() {}
+
+ // These are the virtual functions that must be provided by all
+ // derived patches: GetPou() returns the original smooth
+ // (non-normalized) cutoff on the patch; F() and Inverse() implement
+ // the mapping f: (t)->(x,y,z) and its inverse; and the Df*() and Dn*()
+ // functions return the derivatives of the mapping f and unit normal n
+ // with respect to u and v
+ virtual double GetPou(double t);
+ virtual int GetTag();
+ virtual void F(double t, double &x, double &y, double &z);
+ virtual bool Inverse(double x,double y,double z,double &t);
+};
+
+#endif
diff --git a/contrib/FourierModel/ExactPatch.cpp b/contrib/FourierModel/ExactPatch.cpp
new file mode 100644
index 0000000000..bb69ca7b76
--- /dev/null
+++ b/contrib/FourierModel/ExactPatch.cpp
@@ -0,0 +1,192 @@
+#include "Message.h"
+#include "ExactPatch.h"
+
+ExactPatch::ExactPatch(PatchInfo *PI) : Patch(PI)
+{
+ _NcrossT[0] = _PI->normal[1] * _PI->tangent[2] -
+ _PI->normal[2] * _PI->tangent[1];
+ _NcrossT[1] = _PI->normal[2] * _PI->tangent[0] -
+ _PI->normal[0] * _PI->tangent[2];
+ _NcrossT[2] = _PI->normal[0] * _PI->tangent[1] -
+ _PI->normal[1] * _PI->tangent[0];
+
+ if (!strcmp(_PI->type,"plane")) {
+ _PI->periodic[0] = 0;
+ _PI->periodic[1] = 0;
+ }
+ else if (!strcmp(_PI->type,"disc")) {
+ _PI->periodic[0] = 0;
+ _PI->periodic[1] = 1;
+ }
+ else {
+ Msg::Info("Unknown exact patch type");
+ _PI->periodic[0] = 0;
+ _PI->periodic[1] = 0;
+ }
+}
+
+void ExactPatch::F(double u, double v, double &x, double &y, double &z)
+{
+ if (!strcmp(_PI->type,"plane")) {
+ x = _PI->origin[0] + _PI->tangent[0] * _PI->scale[0] * u +
+ _NcrossT[0] * _PI->scale[1] * v;
+ y = _PI->origin[1] + _PI->tangent[1] * _PI->scale[0] * u +
+ _NcrossT[1] * _PI->scale[1] * v;
+ z = _PI->origin[2] + _PI->tangent[2] * _PI->scale[0] * u +
+ _NcrossT[2] * _PI->scale[1] * v;
+ }
+ else if (!strcmp(_PI->type,"disc")) {
+ x = _PI->origin[0] + _PI->scale[1] * u *
+ (_PI->tangent[0] * cos(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[0] * sin(_PI->scale[0] * (v-0.5)));
+ y = _PI->origin[1] + _PI->scale[1] * u *
+ (_PI->tangent[1] * cos(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[1] * sin(_PI->scale[0] * (v-0.5)));
+ z = _PI->origin[2] + _PI->scale[1] * u *
+ (_PI->tangent[2] * cos(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[2] * sin(_PI->scale[0] * (v-0.5)));
+ }
+ else {
+ Msg::Info("Unknown exact patch type");
+ x = y = z = 0.;
+ }
+}
+
+bool ExactPatch::Inverse(double x,double y,double z,double &u,double &v)
+{
+ if (!strcmp(_PI->type,"plane")) {
+ u = ((x - _PI->origin[0]) * _PI->tangent[0] +
+ (y - _PI->origin[1]) * _PI->tangent[1] +
+ (z - _PI->origin[2]) * _PI->tangent[2]) / _PI->scale[0];
+ v = ((x - _PI->origin[0]) * _NcrossT[0] +
+ (y - _PI->origin[1]) * _NcrossT[1] +
+ (z - _PI->origin[2]) * _NcrossT[2]) / _PI->scale[1];
+ }
+ else if (!strcmp(_PI->type,"disc")) {
+ double a = (x - _PI->origin[0]) * _PI->tangent[0] +
+ (y - _PI->origin[1]) * _PI->tangent[1] +
+ (z - _PI->origin[2]) * _PI->tangent[2];
+ double b = (x - _PI->origin[0]) * _NcrossT[0] +
+ (y - _PI->origin[1]) * _NcrossT[1] + (z - _PI->origin[2]) * _NcrossT[2];
+
+ u = sqrt(a * a + b * b);
+ u /= _PI->scale[1];
+
+ v = atan2(b,a);
+ v /= _PI->scale[0];
+ v += 0.5;
+ }
+ else {
+ Msg::Info("Unknown exat patch type");
+ x = y = z = 0.;
+ }
+ return true;
+}
+
+void ExactPatch::Dfdu(double u, double v, double &x, double &y, double &z)
+{
+ if (!strcmp(_PI->type,"plane")) {
+ x = _PI->tangent[0] * _PI->scale[0];
+ y = _PI->tangent[1] * _PI->scale[0];
+ z = _PI->tangent[2] * _PI->scale[0];
+ }
+ else if (!strcmp(_PI->type,"disc")) {
+ x = _PI->scale[1] * (_PI->tangent[0] * cos(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[0] * sin(_PI->scale[0] * (v-0.5)));
+ y = _PI->scale[1] * (_PI->tangent[1] * cos(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[1] * sin(_PI->scale[0] * (v-0.5)));
+ z = _PI->scale[1] * (_PI->tangent[2] * cos(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[2] * sin(_PI->scale[0] * (v-0.5)));
+ }
+ else {
+ Msg::Info("Unknown exact patch type");
+ x = y = z = 0.;
+ }
+}
+
+void ExactPatch::Dfdv(double u, double v, double &x, double &y, double &z)
+{
+ if (!strcmp(_PI->type,"plane")) {
+ x = _NcrossT[0] * _PI->scale[1];
+ y = _NcrossT[1] * _PI->scale[1];
+ z = _NcrossT[2] * _PI->scale[1];
+ }
+ else if (!strcmp(_PI->type,"disc")) {
+ x = _PI->scale[1] * u * _PI->scale[0] *
+ (- _PI->tangent[0] * sin(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[0] * cos(_PI->scale[0] * (v-0.5)));
+ y = _PI->scale[1] * u * _PI->scale[0] *
+ (- _PI->tangent[1] * sin(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[1] * cos(_PI->scale[0] * (v-0.5)));
+ z = _PI->scale[1] * u * _PI->scale[0] *
+ (- _PI->tangent[2] * sin(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[2] * cos(_PI->scale[0] * (v-0.5)));
+ }
+ else {
+ Msg::Info("Unknown exact patch type");
+ x = y = z = 0.;
+ }
+}
+
+void ExactPatch::Dfdfdudu(double u,double v,double &x,double &y,double &z)
+{
+ if (!strcmp(_PI->type,"plane")) {
+ x = y = z = 0.;
+ }
+ else if (!strcmp(_PI->type,"disc")) {
+ x = y = z = 0.;
+ }
+ else {
+ Msg::Info("Unknown exact patch type");
+ x = y = z = 0.;
+ }
+}
+
+void ExactPatch::Dfdfdudv(double u,double v,double &x,double &y,double &z)
+{
+ if (!strcmp(_PI->type,"plane")) {
+ x = y = z = 0.;
+ }
+ else if (!strcmp(_PI->type,"disc")) {
+ x = _PI->scale[1] * _PI->scale[0] *
+ (- _PI->tangent[0] * sin(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[0] * cos(_PI->scale[0] * (v-0.5)));
+ y = _PI->scale[1] * _PI->scale[0] *
+ (- _PI->tangent[1] * sin(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[1] * cos(_PI->scale[0] * (v-0.5)));
+ z = _PI->scale[1] * _PI->scale[0] *
+ (- _PI->tangent[2] * sin(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[2] * cos(_PI->scale[0] * (v-0.5)));
+ }
+ else {
+ Msg::Info("Unknown exact patch type");
+ x = y = z = 0.;
+ }
+}
+
+void ExactPatch::Dfdfdvdv(double u,double v,double &x,double &y,double &z)
+{
+ if (!strcmp(_PI->type,"plane")) {
+ x = y = z = 0.;
+ }
+ else if (!strcmp(_PI->type,"disc")) {
+ x = - _PI->scale[1] * u * _PI->scale[0] * _PI->scale[0] *
+ (_PI->tangent[0] * cos(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[0] * sin(_PI->scale[0] * (v-0.5)));
+ y = - _PI->scale[1] * u * _PI->scale[0] * _PI->scale[0] *
+ (- _PI->tangent[1] * cos(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[1] * sin(_PI->scale[0] * (v-0.5)));
+ z = - _PI->scale[1] * u * _PI->scale[0] * _PI->scale[0] *
+ (- _PI->tangent[2] * cos(_PI->scale[0] * (v-0.5)) +
+ _NcrossT[2] * sin(_PI->scale[0] * (v-0.5)));
+ }
+ else {
+ Msg::Info("Unknown exact patch type");
+ x = y = z = 0.;
+ }
+}
+
+double ExactPatch::GetPou(double u, double v)
+{
+ return 1.;
+}
diff --git a/contrib/FourierModel/ExactPatch.h b/contrib/FourierModel/ExactPatch.h
new file mode 100644
index 0000000000..dd14469414
--- /dev/null
+++ b/contrib/FourierModel/ExactPatch.h
@@ -0,0 +1,30 @@
+#ifndef _EXACT_PATCH_H_
+#define _EXACT_PATCH_H_
+
+#include "Patch.h"
+
+// The base class for the patches
+class ExactPatch : public Patch {
+ private:
+ // third axis
+ double _NcrossT[3];
+ public:
+ ExactPatch(PatchInfo* PI);
+ virtual ~ExactPatch() {}
+
+ // These are the virtual functions that must be provided by all
+ // derived patches: GetPou() returns the original smooth
+ // (non-normalized) cutoff on the patch; F() and Inverse() implement
+ // the mapping f: (u,v)->(x,y,z) and its inverse; and the Df*() and Dn*()
+ // functions return the derivatives of the mapping f and unit normal n
+ // with respect to u and v
+ virtual double GetPou(double u, double v);
+ virtual void F(double u, double v, double &x, double &y, double &z);
+ virtual bool Inverse(double x,double y,double z,double &u,double &v);
+ virtual void Dfdu(double u, double v, double &x, double &y, double &z);
+ virtual void Dfdv(double u, double v, double &x, double &y, double &z);
+ virtual void Dfdfdudu(double u,double v,double &x,double &y,double &z);
+ virtual void Dfdfdudv(double u,double v,double &x,double &y,double &z);
+ virtual void Dfdfdvdv(double u,double v,double &x,double &y,double &z);
+};
+#endif
diff --git a/contrib/FourierModel/FM_Edge.cpp b/contrib/FourierModel/FM_Edge.cpp
new file mode 100644
index 0000000000..a49f7d5fad
--- /dev/null
+++ b/contrib/FourierModel/FM_Edge.cpp
@@ -0,0 +1,123 @@
+#include "FM_Edge.h"
+#include "Message.h"
+
+void FM_Edge::F(double t, double &x, double &y, double &z)
+{
+ if (_curve) {
+ double tStart, tEnd;
+ _curve->Inverse(_SP->GetX(),_SP->GetY(),_SP->GetZ(),tStart);
+ _curve->Inverse(_EP->GetX(),_EP->GetY(),_EP->GetZ(),tEnd);
+
+ double tRescaled = tStart + t * (tEnd - tStart);
+ _curve->F(tRescaled, x, y, z);
+ //Msg::Info("%g %g %g",_SP->GetX(),_SP->GetY(),_SP->GetZ());
+ //Msg::Info("%g %g %g",_EP->GetX(),_EP->GetY(),_EP->GetZ());
+ //Msg::Info("t : %g %g %g %g",t,tRescaled, tStart, tEnd);
+ }
+ else {
+ x = _SP->GetX() + t * (_EP->GetX() - _SP->GetX());
+ y = _SP->GetY() + t * (_EP->GetY() - _SP->GetY());
+ z = _SP->GetZ() + t * (_EP->GetZ() - _SP->GetZ());
+ }
+}
+
+bool FM_Edge::Inverse(double x,double y,double z,double &t)
+{
+ if (_curve) {
+ double tStart, tEnd;
+ _curve->Inverse(_SP->GetX(),_SP->GetY(),_SP->GetZ(),tStart);
+ _curve->Inverse(_EP->GetX(),_EP->GetY(),_EP->GetZ(),tEnd);
+
+ double tCurve;
+ _curve->Inverse(x, y, z, tCurve);
+ t = (tCurve - tStart) / (tEnd - tStart);
+ }
+ else {
+ if (_EP->GetX() - _SP->GetX())
+ t = (x - _SP->GetX()) / (_EP->GetX() - _SP->GetX());
+ else if (_EP->GetY() - _SP->GetY())
+ t = (y - _SP->GetY()) / (_EP->GetY() - _SP->GetY());
+ else if (_EP->GetZ() - _SP->GetZ())
+ t = (z - _SP->GetZ()) / (_EP->GetZ() - _SP->GetZ());
+ else {
+ Msg::Warning("Cannnot inver the curve");
+ t = 0.;
+ }
+ }
+}
+
+void FM_Edge::Dfdt(double t, double &x, double &y, double &z)
+{
+ if (_curve) {
+ double xMinus, yMinus, zMinus;
+ double xPlus, yPlus, zPlus;
+
+ double tStart, tEnd;
+ _curve->Inverse(_SP->GetX(),_SP->GetY(),_SP->GetZ(),tStart);
+ _curve->Inverse(_EP->GetX(),_EP->GetY(),_EP->GetZ(),tEnd);
+
+ double h = 1.e-10;
+ double tRescaled = tStart + t * (tEnd - tStart);
+ if (t+0.5*h > 1.) {
+ _curve->F(tRescaled, xPlus, yPlus, zPlus);
+ double tMinus = tStart + (t - h) * (tEnd - tStart);
+ _curve->F(tMinus, xMinus, yMinus, zMinus);
+ }
+ else if (t-0.5*h < 0.) {
+ _curve->F(tRescaled, xMinus, yMinus, zMinus);
+ double tPlus = tStart + (t + h) * (tEnd - tStart);
+ _curve->F(tPlus, xPlus, yPlus, zPlus);
+ }
+ else {
+ double tPlus = tStart + (t + 0.5 * h) * (tEnd - tStart);
+ double tMinus = tStart + (t - 0.5 * h) * (tEnd - tStart);
+ _curve->F(tPlus, xPlus, yPlus, zPlus);
+ _curve->F(tMinus, xMinus, yMinus, zMinus);
+ }
+ x = (xPlus - xMinus) / h;
+ y = (yPlus - yMinus) / h;
+ z = (zPlus - zMinus) / h;
+ }
+ else {
+ x = _EP->GetX() - _SP->GetX();
+ y = _EP->GetY() - _SP->GetY();
+ z = _EP->GetZ() - _SP->GetZ();
+ }
+}
+
+void FM_Edge::Dfdfdtdt(double t, double &x, double &y, double &z)
+{
+ if (_curve) {
+ double xMinus, yMinus, zMinus;
+ double xPlus, yPlus, zPlus;
+
+ double tStart, tEnd;
+ _curve->Inverse(_SP->GetX(),_SP->GetY(),_SP->GetZ(),tStart);
+ _curve->Inverse(_EP->GetX(),_EP->GetY(),_EP->GetZ(),tEnd);
+
+ double h = 1.e-10;
+ double tRescaled = tStart + t * (tEnd - tStart);
+ if (t+0.5*h > 1.) {
+ Dfdt(tRescaled, xPlus, yPlus, zPlus);
+ double tMinus = tStart + (t - h) * (tEnd - tStart);
+ Dfdt(tMinus, xMinus, yMinus, zMinus);
+ }
+ else if (t-0.5*h < 0.) {
+ Dfdt(tRescaled, xMinus, yMinus, zMinus);
+ double tPlus = tStart + (t + h) * (tEnd - tStart);
+ Dfdt(tPlus, xPlus, yPlus, zPlus);
+ }
+ else {
+ double tPlus = tStart + (t + 0.5 * h) * (tEnd - tStart);
+ double tMinus = tStart + (t - 0.5 * h) * (tEnd - tStart);
+ Dfdt(tPlus, xPlus, yPlus, zPlus);
+ Dfdt(tMinus, xMinus, yMinus, zMinus);
+ }
+ x = (xPlus - xMinus) / h;
+ y = (yPlus - yMinus) / h;
+ z = (zPlus - zMinus) / h;
+ }
+ else {
+ x = y = z = 0.;
+ }
+}
diff --git a/contrib/FourierModel/FM_Edge.h b/contrib/FourierModel/FM_Edge.h
new file mode 100644
index 0000000000..c80c7b872a
--- /dev/null
+++ b/contrib/FourierModel/FM_Edge.h
@@ -0,0 +1,37 @@
+#ifndef _FM_EDGE_H_
+#define _FM_EDGE_H_
+
+#include "Curve.h"
+#include "FM_Vertex.h"
+
+class FM_Edge {
+ private:
+ Curve* _curve;
+ FM_Vertex* _SP;
+ FM_Vertex* _EP;
+ public:
+ FM_Edge() : _curve(0), _SP(0), _EP(0) {}
+ FM_Edge(Curve* curve) : _curve(curve), _SP(0), _EP(0) {}
+ FM_Edge(Curve* curve, FM_Vertex* SP, FM_Vertex* EP) :
+ _curve(curve), _SP(SP), _EP(EP) {}
+ virtual ~FM_Edge() {}
+
+ inline FM_Vertex* GetStartPoint() { return _SP; }
+ inline FM_Vertex* GetEndPoint() { return _EP; }
+
+ inline void SetStartPoint(FM_Vertex* SP) { _SP = SP; }
+ inline void SetStartPoint(double x, double y, double z) {
+ _SP = new FM_Vertex(x,y,z);
+ }
+ inline void SetEndPoint(FM_Vertex* EP) { _EP = EP; }
+ inline void SetEndPoint(double x, double y, double z) {
+ _EP = new FM_Vertex(x,y,z);
+ }
+
+ void F(double t, double &x, double &y, double &z);
+ bool Inverse(double x,double y,double z,double &t);
+ void Dfdt(double t, double &x, double &y, double &z);
+ void Dfdfdtdt(double t, double &x, double &y, double &z);
+};
+
+#endif
diff --git a/contrib/FourierModel/FM_Face.cpp b/contrib/FourierModel/FM_Face.cpp
new file mode 100644
index 0000000000..a25f26ee02
--- /dev/null
+++ b/contrib/FourierModel/FM_Face.cpp
@@ -0,0 +1,147 @@
+#include "FM_Face.h"
+#include "Message.h"
+
+void FM_Face::F(double u, double v, double &x, double &y, double &z) {
+ if (_edge.size() == 0) {
+ _patch->F(u,v,x,y,z);
+ }
+ else if (_edge.size() == 1) {
+ // sanity check
+ double cx, cy, cz;
+ _patch->F(0.,0.,cx,cy,cz);
+ double px, py, pz;
+ _edge[0]->F(v,px,py,pz);
+ double R = sqrt((px-cx)*(px-cx)+(py-cy)*(py-cy)+(pz-cz)*(pz-cz));
+ _patch->F(u*R,v,x,y,z);
+ }
+ else if (_edge.size() == 4) {
+ double xD, yD, zD, uD, vD;
+ double xU, yU, zU, uU, vU;
+ _edge[0]->F(v,xD,yD,zD);
+ _patch->Inverse(xD,yD,zD,uD,vD);
+ _edge[2]->F(v,xU,yU,zU);
+ _patch->Inverse(xU,yU,zU,uU,vU);
+
+ double U = uD + u * (uU - uD);
+ double V = vD;
+
+ _patch->F(U,V,x,y,z);
+ }
+ else {
+ Msg::Info("Such a face not implemented yet");
+ x = y = z = 0.;
+ }
+}
+
+bool FM_Face::Inverse(double x,double y,double z,double &u, double &v)
+{
+ if (_edge.size() == 0) {
+ _patch->Inverse(x,y,z,u,v);
+ }
+ else if (_edge.size() == 1) {
+ double n[3], t[3], s[3], c[3], r[3];
+ _patch->F(0.,0.,c[0],c[1],c[2]);
+ for (int i=0;i<3;i++) {
+ t[i] = _patch->_PI->tangent[i];
+ n[i] = _patch->_PI->normal[i];
+ }
+
+ s[0] = n[1] * t[2] - n[2] * t[1];
+ s[1] = n[2] * t[0] - n[0] * t[2];
+ s[2] = n[0] * t[1] - n[1] * t[0];
+
+ r[0] = x - c[0];
+ r[1] = y - c[1];
+ r[2] = z - c[2];
+
+ double norm = sqrt(r[0]*r[0]+r[1]*r[1]+r[2]*r[2]);
+ for (int i=0;i<3;i++)
+ r[i] /= norm;
+ double xx, yy;
+ for (int i=0;i<3;i++) {
+ xx += r[i] * t[i];
+ yy += r[i] * s[i];
+ }
+ v = atan2(yy, xx)/(2. * M_PI) +0.5;
+ double px, py, pz;
+ _edge[0]->F(v,px,py,pz);
+ double R = sqrt((px-c[0])*(px-c[0])+(py-c[1])*(py-c[1])+
+ (pz-c[2])*(pz-c[2]));
+ u = norm / R;
+ }
+ else if (_edge.size() == 4) {
+ double U, V;
+ _patch->Inverse(x,y,z,U,V);
+ double u0,v0,u1,v1;
+ double xD, yD, zD, uD, vD;
+ double xU, yU, zU, uU, vU;
+ _edge[0]->F(0,xD,yD,zD);
+ _patch->Inverse(xD,yD,zD,u0,v0);
+ _edge[0]->F(1,xD,yD,zD);
+ _patch->Inverse(xD,yD,zD,u1,v1);
+ v = (V-v0)/(v1-v0);
+ _edge[0]->F(v,xD,yD,zD);
+ _patch->Inverse(xD,yD,zD,uD,vD);
+ _edge[2]->F(v,xU,yU,zU);
+ _patch->Inverse(xU,yU,zU,uU,vU);
+
+ u = (U - uD) / (uU - uD);
+ }
+ else {
+ Msg::Info("Such a face not implemented yet");
+ u = v = 0.;
+ }
+
+}
+
+void FM_Face::Dfdu(double u, double v, double &x, double &y, double &z)
+{
+ Msg::Info("Not implemented yet.");
+ x = y = z = 0.;
+}
+
+void FM_Face::Dfdv(double u, double v, double &x, double &y, double &z)
+{
+ Msg::Info("Not implemented yet.");
+ x = y = z = 0.;
+}
+
+void FM_Face::Dfdfdudu(double u, double v, double &x, double &y, double &z)
+{
+ Msg::Info("Not implemented yet.");
+ x = y = z = 0.;
+}
+
+void FM_Face::Dfdfdudv(double u, double v, double &x, double &y, double &z)
+{
+ Msg::Info("Not implemented yet.");
+ x = y = z = 0.;
+}
+
+void FM_Face::Dfdfdvdv(double u, double v, double &x, double &y, double &z)
+{
+ Msg::Info("Not implemented yet.");
+ x = y = z = 0.;
+}
+
+void FM_Face::GetNormal(double u, double v, double &x, double &y, double &z)
+{
+ double dfdu[3], dfdv[3];
+ Dfdu(u, v, dfdu[0], dfdu[1], dfdu[2]);
+ Dfdv(u, v, dfdv[0], dfdv[1], dfdv[2]);
+ x = dfdu[1] * dfdv[2] - dfdu[2] * dfdv[1];
+ y = dfdu[2] * dfdv[0] - dfdu[0] * dfdv[2];
+ z = dfdu[0] * dfdv[1] - dfdu[1] * dfdv[0];
+}
+
+void FM_Face::
+GetUnitNormal(double u, double v, double &x, double &y, double &z)
+{
+ GetNormal(u, v, x, y, z);
+ double norm = sqrt(x * x + y * y + z * z);
+ if(norm > 0.) {
+ x /= norm;
+ y /= norm;
+ z /= norm;
+ }
+}
diff --git a/contrib/FourierModel/FM_Face.h b/contrib/FourierModel/FM_Face.h
new file mode 100644
index 0000000000..f8f7a6e0c1
--- /dev/null
+++ b/contrib/FourierModel/FM_Face.h
@@ -0,0 +1,33 @@
+#ifndef _FM_FACE_H_
+#define _FM_FACE_H_
+
+#include "Patch.h"
+#include "FM_Edge.h"
+
+class FM_Face {
+ private:
+ Patch* _patch;
+ std::vector<FM_Edge*> _edge;
+ public:
+ FM_Face() : _patch(0) {}
+ FM_Face(Patch* patch) : _patch(patch) {}
+ FM_Face(Patch* patch, std::vector<FM_Edge*> edge) :
+ _patch(patch), _edge(edge) {}
+ virtual ~FM_Face() {}
+
+ inline void AddEdge(FM_Edge* edge) { _edge.push_back(edge); }
+ inline int GetNumEdges() { return _edge.size(); }
+ inline FM_Edge* GetEdge(int i) { return _edge[i]; }
+
+ void F(double u, double v, double &x, double &y, double &z);
+ bool Inverse(double x,double y,double z,double &u, double &v);
+ void Dfdu(double u, double v, double &x, double &y, double &z);
+ void Dfdv(double u, double v, double &x, double &y, double &z);
+ void Dfdfdudu(double u, double v, double &x, double &y, double &z);
+ void Dfdfdudv(double u, double v, double &x, double &y, double &z);
+ void Dfdfdvdv(double u, double v, double &x, double &y, double &z);
+ void GetUnitNormal(double u,double v,double &x,double &y,double &z);
+ void GetNormal(double u, double v, double &x, double &y, double &z);
+};
+
+#endif
diff --git a/contrib/FourierModel/FM_Info.cpp b/contrib/FourierModel/FM_Info.cpp
new file mode 100644
index 0000000000..8b13789179
--- /dev/null
+++ b/contrib/FourierModel/FM_Info.cpp
@@ -0,0 +1 @@
+
diff --git a/contrib/FourierModel/FM_Info.h b/contrib/FourierModel/FM_Info.h
new file mode 100644
index 0000000000..4eb3a0071f
--- /dev/null
+++ b/contrib/FourierModel/FM_Info.h
@@ -0,0 +1,37 @@
+#ifndef _FM_INFO_H_
+#define _FM_INFO_H_
+
+#include <complex>
+#include <vector>
+
+class PatchInfo {
+ public:
+ char type[16];
+ char projection[16];
+ int tag;
+ int nModes[2];
+ int periodic[2];
+ double normal[3];
+ double origin[3];
+ double tangent[3];
+ double scale[2];
+ std::vector<std::vector<std::complex<double> > > coeff;
+
+ PatchInfo() {}
+ virtual ~PatchInfo() {}
+};
+
+class IntersectionInfo {
+ public:
+ int tag;
+ double SP[3];
+ double EP[3];
+ struct {
+ int patchTag;
+ } intersectingPatches[2];
+
+ IntersectionInfo() {}
+ virtual ~IntersectionInfo() {}
+};
+
+#endif
diff --git a/contrib/FourierModel/FM_Reader.cpp b/contrib/FourierModel/FM_Reader.cpp
new file mode 100644
index 0000000000..37cfa6c072
--- /dev/null
+++ b/contrib/FourierModel/FM_Reader.cpp
@@ -0,0 +1,108 @@
+#include <iostream>
+#include <fstream>
+#include "Message.h"
+#include "FM_Reader.h"
+
+FM_Reader::FM_Reader(const char* fn)
+{
+ char c;
+ char continuation[16] = "continuation";
+ std::ifstream InputFile(fn);
+ if (!InputFile) {
+ Msg::Info("Failed to open input file.");
+ exit(EXIT_FAILURE);
+ }
+ InputFile >> _nPatches;
+ InputFile >> _nIntersections;
+
+ for (int i=0;i<_nPatches;i++) {
+ _patchList.push_back(new PatchInfo);;
+ InputFile >> _patchList[i]->tag;
+ InputFile >> _patchList[i]->type;
+ InputFile >> _patchList[i]->projection;
+ InputFile >> _patchList[i]->normal[0] >> _patchList[i]->normal[1] >>
+ _patchList[i]->normal[2];
+ InputFile >> _patchList[i]->origin[0] >> _patchList[i]->origin[1] >>
+ _patchList[i]->origin[2];
+ InputFile >> _patchList[i]->tangent[0] >> _patchList[i]->tangent[1] >>
+ _patchList[i]->tangent[2];
+ InputFile >> _patchList[i]->scale[0] >> _patchList[i]->scale[1];
+ if (!strcmp(_patchList[i]->type,continuation)) {
+ InputFile >> _patchList[i]->periodic[0] >> _patchList[i]->periodic[1];
+ InputFile >> _patchList[i]->nModes[0] >> _patchList[i]->nModes[1];
+ _patchList[i]->coeff.resize(_patchList[i]->nModes[0]);
+ for (int j=0;j<_patchList[i]->nModes[0];j++) {
+ _patchList[i]->coeff[j].resize(_patchList[i]->nModes[1]);
+ for (int k=0;k<_patchList[i]->nModes[1];k++) {
+ double realCoeff, imagCoeff;
+ InputFile >> realCoeff >> imagCoeff;
+ _patchList[i]->coeff[j][k] =
+ std::complex<double>(realCoeff,imagCoeff);
+ }
+ }
+ }
+ }
+ for (int i=0;i<_nIntersections;i++) {
+ _intersectionList.push_back(new IntersectionInfo);
+ InputFile >> _intersectionList[i]->tag;
+ InputFile >> _intersectionList[i]->SP[0] >> _intersectionList[i]->SP[1]
+ >> _intersectionList[i]->SP[2];
+ InputFile >> _intersectionList[i]->EP[0] >> _intersectionList[i]->EP[1]
+ >> _intersectionList[i]->EP[2];
+ InputFile >> _intersectionList[i]->intersectingPatches[0].patchTag;
+ InputFile >> _intersectionList[i]->intersectingPatches[1].patchTag;
+ }
+
+ _patch.resize(_nPatches);
+ _intersection.resize(_nIntersections);
+
+ for (int i=0;i<_nPatches;i++) {
+ PatchInfo* PI = _patchList[i];
+ if (!strcmp(PI->type,continuation))
+ _patch[PI->tag] = new ContinuationPatch(PI,2);
+ else
+ _patch[PI->tag] = new ExactPatch(PI);
+ }
+ for (int i=0;i<_nIntersections;i++) {
+ IntersectionInfo* II = _intersectionList[i];
+ _intersection[II->tag] = new Curve(II,_patch);
+ }
+
+ InputFile >> _nFaces;
+ for (int i=0;i<_nFaces;i++) {
+ int faceTag, nEdges;
+ InputFile >> faceTag;
+ _face.push_back(new FM_Face(GetPatch(faceTag)));
+ InputFile >> nEdges;
+ for (int j=0;j<nEdges;j++) {
+ int edgeTag;
+ double SPx,SPy,SPz;
+ double EPx,EPy,EPz;
+ InputFile >> edgeTag;
+ InputFile >> SPx >> SPy >> SPz;
+ InputFile >> EPx >> EPy >> EPz;
+ if (edgeTag < 0)
+ _face[i]->AddEdge(new FM_Edge(0,
+ new FM_Vertex(SPx,SPy,SPz),
+ new FM_Vertex(EPx,EPy,EPz)));
+ else
+ _face[i]->AddEdge(new FM_Edge(GetIntersection(edgeTag),
+ new FM_Vertex(SPx,SPy,SPz),
+ new FM_Vertex(EPx,EPy,EPz)));
+ }
+ }
+}
+
+Patch* FM_Reader::GetPatch(int tag)
+{
+ for (int i=0;i<_patch.size();i++)
+ if (_patch[i]->GetTag() == tag)
+ return _patch[i];
+}
+
+Curve* FM_Reader::GetIntersection(int tag)
+{
+ for (int i=0;i<_intersection.size();i++)
+ if (_intersection[i]->GetTag() == tag)
+ return _intersection[i];
+}
diff --git a/contrib/FourierModel/FM_Reader.h b/contrib/FourierModel/FM_Reader.h
new file mode 100644
index 0000000000..98940271bf
--- /dev/null
+++ b/contrib/FourierModel/FM_Reader.h
@@ -0,0 +1,38 @@
+#ifndef _FM_READER_H_
+#define _FM_READER_H_
+
+#include <vector>
+#include <complex>
+#include "Curve.h"
+#include "ExactPatch.h"
+#include "ContinuationPatch.h"
+#include "FM_Info.h"
+#include "FM_Face.h"
+
+class FM_Reader {
+ private:
+ int _nFaces;
+ int _nPatches;
+ int _nIntersections;
+ std::vector<PatchInfo*> _patchList;
+ std::vector<IntersectionInfo*> _intersectionList;
+ std::vector<FM_Face*> _face;
+ std::vector<Patch*> _patch;
+ std::vector<Curve*> _intersection;
+ public:
+ FM_Reader(const char* fn);
+ virtual ~FM_Reader() {}
+
+ inline int GetNumFaces() { return _nFaces; }
+ inline int GetNumPatches() { return _nPatches; }
+ inline int GetNumIntersections() { return _nIntersections; }
+ inline std::vector<PatchInfo*> GetPatchList() { return _patchList; }
+ inline std::vector<IntersectionInfo*> GetIntersectionList()
+ { return _intersectionList; }
+ inline FM_Face* GetFace(int i) { return _face[i]; }
+
+ Patch* GetPatch(int tag);
+ Curve* GetIntersection(int tag);
+};
+
+#endif
diff --git a/contrib/FourierModel/FM_Vertex.cpp b/contrib/FourierModel/FM_Vertex.cpp
new file mode 100644
index 0000000000..139597f9cb
--- /dev/null
+++ b/contrib/FourierModel/FM_Vertex.cpp
@@ -0,0 +1,2 @@
+
+
diff --git a/contrib/FourierModel/FM_Vertex.h b/contrib/FourierModel/FM_Vertex.h
new file mode 100644
index 0000000000..2bc0fe0cae
--- /dev/null
+++ b/contrib/FourierModel/FM_Vertex.h
@@ -0,0 +1,21 @@
+#ifndef _FM_VERTEX_H_
+#define _FM_VERTEX_H_
+
+class FM_Vertex {
+ private:
+ double _x,_y,_z;
+ public:
+ FM_Vertex() : _x(0), _y(0), _z(0) {}
+ FM_Vertex(double x, double y, double z) : _x(x), _y(y), _z(z) {}
+ virtual ~FM_Vertex() {}
+
+ inline double GetX() { return _x; }
+ inline double GetY() { return _y; }
+ inline double GetZ() { return _z; }
+
+ inline void SetX(double x) { _x = x; }
+ inline void SetY(double y) { _y = y; }
+ inline void SetZ(double z) { _z = z; }
+};
+
+#endif
diff --git a/contrib/FourierModel/Main.cpp b/contrib/FourierModel/Main.cpp
new file mode 100644
index 0000000000..1dc9e6185a
--- /dev/null
+++ b/contrib/FourierModel/Main.cpp
@@ -0,0 +1,117 @@
+#include <cstring>
+#include <iostream>
+#include "Utils.h"
+#include "Message.h"
+#include "FM_Reader.h"
+
+int main(int argc, char *argv[])
+{
+ char* fn;
+ char continuation[16] = "continuation";
+
+ if (argc == 1) {
+ Msg::Info("Reading from default file : untitled.fm");
+ fn = "untitled.fm";
+ }
+ else
+ fn = argv[1];
+
+ FM_Reader reader(fn);
+
+ /*
+
+ int n1=64;
+ int n2=64;
+ double h1 = 1./(n1-1);
+ double h2 = 1./(n2-1);
+ for (int j=0;j<n1;j++) {
+ for (int k=0;k<n2;k++) {
+ double u = j*h1;
+ double v = k*h2;
+ double x, y, z;
+ patches[0]->F(u,v,x,y,z);
+ printf("x(%d,%d) = %g; y(%d,%d) = %g; z(%d,%d) = %g;\n",
+ j+1,k+1,x,j+1,k+1,y,j+1,k+1,z);
+ }
+ }
+
+ for (int i=0;i<reader.GetNumIntersections();i++) {
+ int n=128;
+ double h = 1./(n-1);
+ for (int j=0;j<n;j++) {
+ double t = j*h;
+ double x, y, z;
+ reader.GetIntersection(i)->F(t,x,y,z);
+ printf("x(%d,%d) = %g; y(%d,%d) = %g; z(%d,%d) = %g;\n",
+ i+1,j+1,x,i+1,j+1,y,i+1,j+1,z);
+ }
+ }
+
+ for (int i=0;i<reader.GetNumPatches();i++) {
+ std::cout << reader.GetPatchList()[i]->tag << "\n";
+ std::cout << reader.GetPatchList()[i]->type << "\n";
+ std::cout << reader.GetPatchList()[i]->projection << "\n";
+ std::cout << reader.GetPatchList()[i]->normal[0] << " " <<
+ reader.GetPatchList()[i]->normal[1] << " " <<
+ reader.GetPatchList()[i]->normal[2] << "\n";
+ std::cout << reader.GetPatchList()[i]->origin[0] << " " <<
+ reader.GetPatchList()[i]->origin[1] << " " <<
+ reader.GetPatchList()[i]->origin[2] << "\n";
+ std::cout << reader.GetPatchList()[i]->tangent[0] << " " <<
+ reader.GetPatchList()[i]->tangent[1] << " " <<
+ reader.GetPatchList()[i]->tangent[2] << "\n";
+ if (!strcmp(reader.GetPatchList()[i]->type,continuation)) {
+ std::cout << reader.GetPatchList()[i]->periodic[0] << " " <<
+ reader.GetPatchList()[i]->periodic[1] << "\n";
+ std::cout << reader.GetPatchList()[i]->nModes[0] << " " <<
+ reader.GetPatchList()[i]->nModes[1] << "\n";
+ for (int j=0;j<reader.GetPatchList()[i]->nModes[0];j++)
+ for (int k=0;k<reader.GetPatchList()[i]->nModes[1];k++)
+ std::cout << reader.GetPatchList()[i]->coeff[j][k].real() << " " <<
+ reader.GetPatchList()[i]->coeff[j][k].imag() << "\n";
+ }
+ }
+
+ for (int i=0;i<reader.GetIntersectionList().size();i++) {
+ std::cout << reader.GetIntersectionList()[i]->tag << "\n";
+ std::cout << reader.GetIntersectionList()[i]->SP[0] <<
+ " " << reader.GetIntersectionList()[i]->SP[1] << " " <<
+ reader.GetIntersectionList()[i]->SP[2] << "\n";
+ std::cout << reader.GetIntersectionList()[i]->EP[0] <<
+ " " << reader.GetIntersectionList()[i]->EP[1] << " " <<
+ reader.GetIntersectionList()[i]->EP[2] << "\n";
+ std::cout << reader.GetIntersectionList()[i]->
+ intersectingPatches[0].patchTag << "\n";
+ std::cout << reader.GetIntersectionList()[i]->
+ intersectingPatches[1].patchTag << "\n";
+ }
+ */
+
+ int nU=64;
+ int nV=64;
+
+ std::vector<int> color(3);
+ color[0] = 0; color[1] = 0; color[2] = 1;
+
+ std::vector<std::vector<double> > x(nU,std::vector<double>(nV));
+ std::vector<std::vector<double> > y(nU,std::vector<double>(nV));
+ std::vector<std::vector<double> > z(nU,std::vector<double>(nV));
+
+ std::vector<std::vector<int> > mask = ones(nU,nV);
+
+ for (int i=0;i<reader.GetNumFaces();i++) {
+ double hU = 1./(nU-1);
+ double hV = 1./(nV-1);
+ for (int j=0;j<nU;j++) {
+ for (int k=0;k<nV;k++) {
+ double u = j*hU;
+ double v = k*hV;
+ reader.GetFace(i)->F(u,v,x[j][k],y[j][k],z[j][k]);
+ }
+ }
+ if (i == 0)
+ plotSceneViewer(0,"snc.iv",color,x,y,z,nU,nV,mask);
+ else
+ plotSceneViewer(1,"snc.iv",color,x,y,z,nU,nV,mask);
+ }
+}
diff --git a/contrib/FourierModel/Makefile b/contrib/FourierModel/Makefile
new file mode 100644
index 0000000000..5e33ac5b8c
--- /dev/null
+++ b/contrib/FourierModel/Makefile
@@ -0,0 +1,55 @@
+include ../../variables
+
+LIB = ../../lib/libGmshFourierModel.a
+
+CFLAGS = ${OPTIM} ${FLAGS}
+
+SRC = ContinuationPatch.cpp\
+ ExactPatch.cpp\
+ Patch.cpp\
+ Curve.cpp\
+ FM_Edge.cpp\
+ FM_Face.cpp\
+ FM_Info.cpp\
+ FM_Reader.cpp\
+ FM_Vertex.cpp\
+ Message.cpp
+
+OBJ = ${SRC:.cpp=.o}
+
+.SUFFIXES: .o .cpp
+
+${LIB}: ${OBJ}
+ ${AR} ${LIB} ${OBJ}
+ ${RANLIB} ${LIB}
+
+.cpp.o:
+ ${CXX} ${CFLAGS} -c $<
+
+clean:
+ rm -f *.o
+
+depend:
+ (sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \
+ ${CXX} -MM ${CFLAGS} ${SRC} \
+ ) >Makefile.new
+ cp Makefile Makefile.bak
+ cp Makefile.new Makefile
+ rm -f Makefile.new
+
+# DO NOT DELETE THIS LINE
+ContinuationPatch.o: ContinuationPatch.cpp Message.h ContinuationPatch.h \
+ Patch.h FM_Info.h
+ExactPatch.o: ExactPatch.cpp Message.h ExactPatch.h Patch.h FM_Info.h
+Patch.o: Patch.cpp Patch.h FM_Info.h
+Curve.o: Curve.cpp Message.h Curve.h Patch.h FM_Info.h
+FM_Edge.o: FM_Edge.cpp FM_Edge.h Curve.h Patch.h FM_Info.h FM_Vertex.h \
+ Message.h
+FM_Face.o: FM_Face.cpp FM_Face.h Patch.h FM_Info.h FM_Edge.h Curve.h \
+ FM_Vertex.h Message.h
+FM_Info.o: FM_Info.cpp
+FM_Reader.o: FM_Reader.cpp Message.h FM_Reader.h Curve.h Patch.h \
+ FM_Info.h ExactPatch.h ContinuationPatch.h FM_Face.h FM_Edge.h \
+ FM_Vertex.h
+FM_Vertex.o: FM_Vertex.cpp
+Message.o: Message.cpp Message.h
diff --git a/contrib/FourierModel/Message.cpp b/contrib/FourierModel/Message.cpp
new file mode 100644
index 0000000000..2034750c02
--- /dev/null
+++ b/contrib/FourierModel/Message.cpp
@@ -0,0 +1,156 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <sys/time.h>
+#include <sys/resource.h>
+
+#if defined(__APPLE__)
+#define RUSAGE_SELF 0
+#define RUSAGE_CHILDREN -1
+#endif
+
+#include "Message.h"
+
+int Message::_commRank = 0;
+int Message::_commSize = 1;
+int Message::_verbosity = 3;
+int Message::_progressMeterStep = 10;
+int Message::_progressMeterCurrent = 0;
+
+void Message::Fatal(char *fmt, ...)
+{
+ va_list args;
+ va_start(args, fmt);
+ if(_commSize > 1)
+ fprintf(stderr, "Fatal : [On processor %d] ", _commRank);
+ else
+ fprintf(stderr, "Fatal : ");
+ vfprintf(stderr, fmt, args);
+ fprintf(stderr, "\n");
+ va_end(args);
+ exit(1);
+}
+
+void Message::Error(char *fmt, ...)
+{
+ va_list args;
+ va_start(args, fmt);
+ if(_commSize > 1)
+ fprintf(stderr, "Error : [On processor %d] ", _commRank);
+ else
+ fprintf(stderr, "Error : ");
+ vfprintf(stderr, fmt, args);
+ fprintf(stderr, "\n");
+ va_end(args);
+}
+
+void Message::Warning(char *fmt, ...)
+{
+ if(_commRank) return;
+ if(_verbosity >= 1){
+ va_list args;
+ va_start(args, fmt);
+ fprintf(stderr, "Warning : ");
+ vfprintf(stderr, fmt, args);
+ fprintf(stderr, "\n");
+ va_end(args);
+ }
+}
+
+void Message::Info(char *fmt, ...)
+{
+ if(_commRank) return;
+ if(_verbosity >= 2){
+ va_list args;
+ va_start(args, fmt);
+ fprintf(stderr, "Info : ");
+ vfprintf(stderr, fmt, args);
+ fprintf(stderr, "\n");
+ va_end(args);
+ }
+}
+
+void Message::Debug(char *fmt, ...)
+{
+ if(_commRank) return;
+ if(_verbosity >= 99){
+ va_list args;
+ va_start(args, fmt);
+ fprintf(stderr, "Debug : ");
+ vfprintf(stderr, fmt, args);
+ fprintf(stderr, "\n");
+ va_end(args);
+ }
+}
+
+void Message::Cpu(char *fmt, ...)
+{
+ if(_commRank) return;
+ if(_verbosity >= 1){
+ static struct rusage r;
+ getrusage(RUSAGE_SELF, &r);
+ double s = (double)r.ru_utime.tv_sec;
+ double us = (double)r.ru_utime.tv_usec;
+ double mem = (double)r.ru_maxrss;
+ va_list args;
+ va_start(args, fmt);
+ fprintf(stderr, "Info : ");
+ vfprintf(stderr, fmt, args);
+ if(mem)
+ fprintf(stderr, " (CPU = %gs Mem = %gkb)\n", s + 1.e-6 * us, mem);
+ else
+ fprintf(stderr, " (CPU = %gs)\n", s + 1.e-6 * us);
+ va_end(args);
+ }
+}
+
+void Message::ResetProgressMeter(int step)
+{
+ _progressMeterStep = step;
+ _progressMeterCurrent = 0;
+}
+
+void Message::ProgressMeter(int n, int N, char *fmt, ...)
+{
+ if(_commRank) return;
+ if(_verbosity >= 2){
+ va_list args;
+ va_start(args, fmt);
+ if(100. * (double)n/(double)N >= _progressMeterCurrent){
+ vfprintf(stderr, fmt, args);
+ fprintf(stderr, "(%d %%) \r", _progressMeterCurrent);
+ _progressMeterCurrent += _progressMeterStep;
+ }
+ if(n >= N - 1)
+ fprintf(stderr, "Done! \r");
+ }
+}
+
+void Message::ProgressMeter(int n, int N)
+{
+ if(_commRank) return;
+ if(_verbosity >= 2){
+ if(100. * (double)n/(double)N >= _progressMeterCurrent){
+ fprintf(stderr, "(%d %%)\r", _progressMeterCurrent);
+ _progressMeterCurrent += _progressMeterStep;
+ }
+ if(n >= N - 1)
+ fprintf(stderr, " \r") ;
+ }
+}
+
+#if defined(HAVE_PETSC)
+
+#include "petsc.h"
+
+void Message::Barrier()
+{
+ MPI_Barrier(PETSC_COMM_WORLD);
+}
+
+#else
+
+void Message::Barrier()
+{
+}
+
+#endif
diff --git a/contrib/FourierModel/Message.h b/contrib/FourierModel/Message.h
new file mode 100644
index 0000000000..f5b9a65464
--- /dev/null
+++ b/contrib/FourierModel/Message.h
@@ -0,0 +1,36 @@
+#ifndef _MESSAGE_H_
+#define _MESSAGE_H_
+
+#include <stdarg.h>
+
+// a class to manage messages
+class Message {
+ private:
+ // current cpu number and total number of cpus
+ static int _commRank, _commSize;
+ // verbosity level
+ static int _verbosity;
+ // step (in %) of the progress meter and current progress %
+ static int _progressMeterStep, _progressMeterCurrent;
+ public:
+ Message() {}
+ static int GetCommRank(){ return _commRank; }
+ static int GetCommSize(){ return _commSize; }
+ static void SetCommRank(int val){ _commRank = val; }
+ static void SetCommSize(int val){ _commSize = val; }
+ static void Barrier();
+ static void SetVerbosity(int val){ _verbosity = val; }
+ static void Fatal(char *fmt, ...);
+ static void Error(char *fmt, ...);
+ static void Warning(char *fmt, ...);
+ static void Info(char *fmt, ...);
+ static void Debug(char *fmt, ...);
+ static void Cpu(char *fmt, ...);
+ static void ProgressMeter(int n, int N);
+ static void ProgressMeter(int n, int N, char *fmt, ...);
+ static void ResetProgressMeter(int step=10);
+};
+
+typedef Message Msg;
+
+#endif
diff --git a/contrib/FourierModel/Patch.cpp b/contrib/FourierModel/Patch.cpp
new file mode 100644
index 0000000000..c8bf1ec899
--- /dev/null
+++ b/contrib/FourierModel/Patch.cpp
@@ -0,0 +1,84 @@
+#include <cmath>
+#include "Patch.h"
+
+void Patch::GetNormal(double u, double v, double &x, double &y, double &z)
+{
+ double dfdu[3], dfdv[3];
+ Dfdu(u, v, dfdu[0], dfdu[1], dfdu[2]);
+ Dfdv(u, v, dfdv[0], dfdv[1], dfdv[2]);
+ x = dfdu[1] * dfdv[2] - dfdu[2] * dfdv[1];
+ y = dfdu[2] * dfdv[0] - dfdu[0] * dfdv[2];
+ z = dfdu[0] * dfdv[1] - dfdu[1] * dfdv[0];
+}
+
+void Patch::GetUnitNormal(double u, double v, double &x, double &y, double &z)
+{
+ GetNormal(u, v, x, y, z);
+ double norm = sqrt(x * x + y * y + z * z);
+ if(norm > 0.) {
+ x /= norm;
+ y /= norm;
+ z /= norm;
+ }
+}
+
+void Patch::Dndu(double u, double v, double &x, double &y, double &z)
+{
+ double n[3],dfdu[3],dfdv[3],dfdfdudu[3],dfdfdudv[3],dfdfdvdv[3],dndu[3];
+ Dfdu(u, v, dfdu[0], dfdu[1], dfdu[2]);
+ Dfdv(u, v, dfdv[0], dfdv[1], dfdv[2]);
+ Dfdfdudu(u, v, dfdfdudu[0], dfdfdudu[1], dfdfdudu[2]);
+ Dfdfdudv(u, v, dfdfdudv[0], dfdfdudv[1], dfdfdudv[2]);
+ Dfdfdvdv(u, v, dfdfdvdv[0], dfdfdvdv[1], dfdfdvdv[2]);
+ GetNormal(u,v,n[0],n[1],n[2]);
+ double norm = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
+ double normSquared = norm * norm;
+ double normCubed = normSquared * norm;
+
+ dndu[0] = dfdfdudu[1] * dfdv[2] + dfdu[1] * dfdfdudv[2] -
+ dfdfdudu[2] * dfdv[1] - dfdu[2] * dfdfdudv[1];
+ dndu[1] = dfdfdudu[2] * dfdv[0] + dfdu[2] * dfdfdudv[0] -
+ dfdfdudu[0] * dfdv[2] - dfdu[0] * dfdfdudv[2];
+ dndu[2] = dfdfdudu[0] * dfdv[1] + dfdu[0] * dfdfdudv[1] -
+ dfdfdudu[1] * dfdv[0] - dfdu[1] * dfdfdudv[0];
+
+ x = ((normSquared - n[0]*n[0])*dndu[0] - n[0]*n[1]*dndu[1] -
+ n[0]*n[2]*dndu[2]) / normCubed;
+ y = ((normSquared - n[1]*n[1])*dndu[1] - n[1]*n[0]*dndu[0] -
+ n[1]*n[2]*dndu[2]) / normCubed;
+ z = ((normSquared - n[2]*n[2])*dndu[2] - n[2]*n[0]*dndu[0] -
+ n[2]*n[1]*dndu[1]) / normCubed;
+}
+
+void Patch::Dndv(double u, double v, double &x, double &y, double &z)
+{
+ double n[3],dfdu[3],dfdv[3],dfdfdudu[3],dfdfdudv[3],dfdfdvdv[3],dndv[3];
+ Dfdu(u, v, dfdu[0], dfdu[1], dfdu[2]);
+ Dfdv(u, v, dfdv[0], dfdv[1], dfdv[2]);
+ Dfdfdudu(u, v, dfdfdudu[0], dfdfdudu[1], dfdfdudu[2]);
+ Dfdfdudv(u, v, dfdfdudv[0], dfdfdudv[1], dfdfdudv[2]);
+ Dfdfdvdv(u, v, dfdfdvdv[0], dfdfdvdv[1], dfdfdvdv[2]);
+ GetNormal(u,v,n[0],n[1],n[2]);
+ double norm = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
+ double normSquared = norm * norm;
+ double normCubed = normSquared * norm;
+
+ dndv[0] = dfdfdudv[1] * dfdv[2] + dfdu[1] * dfdfdvdv[2] -
+ dfdfdudv[2] * dfdv[1] - dfdu[2] * dfdfdvdv[1];
+ dndv[1] = dfdfdudv[2] * dfdv[0] + dfdu[2] * dfdfdvdv[0] -
+ dfdfdudv[0] * dfdv[2] - dfdu[0] * dfdfdvdv[2];
+ dndv[2] = dfdfdudv[0] * dfdv[1] + dfdu[0] * dfdfdvdv[1] -
+ dfdfdudv[1] * dfdv[0] - dfdu[1] * dfdfdvdv[0];
+
+ x = ((normSquared - n[0]*n[0])*dndv[0] - n[0]*n[1]*dndv[1] -
+ n[0]*n[2]*dndv[2]) / normCubed;
+ y = ((normSquared - n[1]*n[1])*dndv[1] - n[1]*n[0]*dndv[0] -
+ n[1]*n[2]*dndv[2]) / normCubed;
+ z = ((normSquared - n[2]*n[2])*dndv[2] - n[2]*n[0]*dndv[0] -
+ n[2]*n[1]*dndv[1]) / normCubed;
+}
+
+int Patch::GetTag()
+{
+ return _PI->tag;
+}
diff --git a/contrib/FourierModel/Patch.h b/contrib/FourierModel/Patch.h
new file mode 100644
index 0000000000..60d78a46b0
--- /dev/null
+++ b/contrib/FourierModel/Patch.h
@@ -0,0 +1,40 @@
+#ifndef _PATCH_H_
+#define _PATCH_H_
+
+#include "FM_Info.h"
+
+// The base class for the patches
+class Patch {
+ public:
+ PatchInfo* _PI;
+ Patch() :_PI(0) {}
+ Patch(PatchInfo* PI) : _PI(PI) {}
+ virtual ~Patch() {}
+
+ int GetTag();
+
+ // These are the virtual functions that must be provided by all
+ // derived patches: GetPou() returns the original smooth
+ // (non-normalized) cutoff on the patch; F() and Inverse() implement
+ // the mapping f: (u,v)->(x,y,z) and its inverse; and the Df*() and Dn*()
+ // functions return the derivatives of the mapping f and unit normal n
+ // with respect to u and v
+ virtual double GetPou(double u, double v) = 0;
+ virtual void F(double u, double v, double &x, double &y, double &z) = 0;
+ virtual bool Inverse(double x,double y,double z,double &u,double &v) = 0;
+ virtual void Dfdu(double u, double v, double &x, double &y, double &z) = 0;
+ virtual void Dfdv(double u, double v, double &x, double &y, double &z) = 0;
+ virtual void Dfdfdudu(double u,double v,double &x,double &y,double &z) = 0;
+ virtual void Dfdfdudv(double u,double v,double &x,double &y,double &z) = 0;
+ virtual void Dfdfdvdv(double u,double v,double &x,double &y,double &z) = 0;
+
+ // These functions may also be provided by the derived patches
+ // (usually for better performance), but they don't have to
+ virtual void Dndu(double u, double v, double &x, double &y, double &z);
+ virtual void Dndv(double u, double v, double &x, double &y, double &z);
+ virtual void GetUnitNormal(double u,double v,double &x,double &y,double &z);
+ virtual void GetNormal(double u, double v, double &x, double &y, double &z);
+
+};
+
+#endif
--
GitLab