diff --git a/Geo/FFace.cpp b/Geo/FFace.cpp
index 043191f86dac5cc0e2d4881a2ca1c105fb8a9b1a..5f391f5742962cfa73cf4693d6bd01b69a1827a2 100644
--- a/Geo/FFace.cpp
+++ b/Geo/FFace.cpp
@@ -38,6 +38,23 @@ 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);
+ printf("per : %d %d\n",face->IsPeriodic(0),face->IsPeriodic(1));
+ if (face->IsPeriodic(0)) {
+ double s,e;
+ printf("\tper : %d %d\n",face->GetEdge(0)->GetCurveExtent(s,e),
+ face->GetEdge(1)->GetCurveExtent(s,e));
+ if (face->GetEdge(0)->GetCurveExtent(s,e)) {
+ printf("%g %g %g\n",u,s,e);
+ u -= floor(u - s);
+ printf("%g %g %g\n\n",u,s,e);
+ }
+ }
+ if (face->IsPeriodic(1)) {
+ double s,e;
+ if (face->GetEdge(1)->GetCurveExtent(s,e)) {
+ v -= floor(v - s);
+ }
+ }
return SPoint2(u, v);
}
@@ -71,4 +88,15 @@ GEntity::GeomType FFace::geomType() const
return GEntity::ParametricSurface;
}
+bool FFace::surfPeriodic(int dim) const
+{
+ return face->IsPeriodic(dim);
+}
+
+double FFace::period(int dir) const
+{
+ return 1.;
+}
+
+
#endif
diff --git a/Geo/FFace.h b/Geo/FFace.h
index fcdf396c077f76c5f7aaa2ecc60427feb84b7a8d..452d26142d2c6efdfd45f240ea06be4a4f11ac55 100644
--- a/Geo/FFace.h
+++ b/Geo/FFace.h
@@ -13,7 +13,6 @@
class FFace : public GFace {
protected:
FM_Face *face;
- bool _periodic[2];
public:
FFace(GModel *m, FM_Face *face_, int tag, std::list<GEdge*> l_edges_);
@@ -37,10 +36,10 @@ class FFace : public GFace {
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;}
+ virtual double period(int dir) const;
ModelType getNativeType() const { return FourierModel; }
void * getNativePtr() const {throw;}
- virtual bool surfPeriodic(int dim) const {throw;}
+ virtual bool surfPeriodic(int dim) const;
};
#endif
diff --git a/Geo/FProjectionFace.cpp b/Geo/FProjectionFace.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..84bc943c36a2583020d2caab7d238af89e100228
--- /dev/null
+++ b/Geo/FProjectionFace.cpp
@@ -0,0 +1,67 @@
+#include "FProjectionFace.h"
+
+#if defined(HAVE_FOURIER_MODEL)
+
+FProjectionFace::FProjectionFace(GModel *m, int num)
+ : GFace(m,num), ps_(0) {}
+
+FProjectionFace::
+FProjectionFace(GModel *m, int num, ProjectionSurface* ps)
+ : GFace(m,num), ps_(ps) {}
+
+FProjectionFace::~FProjectionFace() {}
+
+/*
+void FProjectionFace::ResetProjectionSurface(int type)
+{
+ if (ps_->GetTag() != type) {
+ delete ps_;
+ if (type == 1)
+ ps_ = new RevolvedParabolaProjectionSurface(type);
+ else
+ ps_ = new CylindricalProjectionSurface(type);
+ }
+}
+*/
+
+GPoint FProjectionFace::point(double par1, double par2) const
+{
+ SVector3 p;
+
+ ps_->F(par1,par2,p[0],p[1],p[2]);
+
+ return GPoint(p[0],p[1],p[2]);
+}
+
+SPoint2 FProjectionFace::parFromPoint(const SPoint3 &p) const
+{
+ double u,v;
+ ps_->Inverse(p[0],p[1],p[2],u,v);
+ return SPoint2(u,v);
+}
+
+Pair<SVector3,SVector3> FProjectionFace::firstDer
+(const SPoint2 ¶m) const
+{
+ SVector3 du;
+ SVector3 dv;
+
+ ps_->Dfdu(param.x(),param.y(),du[0],du[1],du[2]);
+ ps_->Dfdv(param.x(),param.y(),dv[0],dv[1],dv[2]);
+
+ return Pair<SVector3,SVector3>(du,dv);
+}
+
+SVector3 FProjectionFace::normal(const SPoint2 ¶m) const
+{
+ double x,y,z;
+ ps_->GetUnitNormal(param.x(),param.y(),x,y,z);
+ return SVector3(x,y,z);
+}
+
+Range<double> FProjectionFace::parBounds(int i) const
+{
+ return Range<double>(0, 1);
+}
+
+#endif
diff --git a/Geo/FProjectionFace.h b/Geo/FProjectionFace.h
new file mode 100644
index 0000000000000000000000000000000000000000..2afd96b2750988089563ffb6ee403f2a1fa2ce0d
--- /dev/null
+++ b/Geo/FProjectionFace.h
@@ -0,0 +1,51 @@
+#ifndef _F_PROJECTION_FACE_H_
+#define _F_PROJECTION_FACE_H_
+
+#include "GModel.h"
+#include "Range.h"
+#include "projectionFace.h"
+#include "CylindricalProjectionSurface.h"
+#include "RevolvedParabolaProjectionSurface.h"
+
+#if defined(HAVE_FOURIER_MODEL)
+
+#include "ProjectionSurface.h"
+
+class FProjectionFace : public GFace {
+ protected:
+ ProjectionSurface* ps_;
+
+ public:
+ FProjectionFace(GModel *m, int num);
+ FProjectionFace(GModel *m, int num, ProjectionSurface* ps);
+ ~FProjectionFace();
+
+ Range<double> parBounds(int i) const;
+ GPoint point(double par1, double par2) const;
+ SVector3 normal(const SPoint2 ¶m) const;
+ Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ SPoint2 parFromPoint(const SPoint3 &) const;
+
+ virtual int paramDegeneracies(int dir, double *par) { return 0; }
+ virtual GPoint closestPoint(const SPoint3 & queryPoint) const {throw;}
+ virtual int containsPoint(const SPoint3 &pt) const {throw;}
+ virtual int containsParam(const SPoint2 &pt) const {throw;}
+ virtual GEntity::GeomType geomType() const
+ { return GEntity::ProjectionFace; }
+ 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 UnknownModel; }
+ void * getNativePtr() const {throw;}
+ virtual bool surfPeriodic(int dim) const {throw;}
+
+ inline ProjectionSurface* GetProjectionSurface() { return ps_; }
+
+ //virtual void ResetProjectionSurface(int type);
+};
+
+#endif
+
+#endif
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index ae4dc268764e02d18980f4bae99a7e67b066dd02..0a45583fcceb29af4a9ae5c34d32ac866674253b 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: GEdge.cpp,v 1.29 2007-05-24 14:44:06 remacle Exp $
+// $Id: GEdge.cpp,v 1.30 2007-07-26 16:28:27 anand Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 50ad8d0e8999d961b8f27e4d8816d04e338ae6f6..9eceb2938aaaabfac4a194013f2eeffeafa309c3 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -82,7 +82,7 @@ class GEntity {
Torus,
RuledSurface,
ParametricSurface,
- ProjectionSurface,
+ ProjectionFace,
BSplineSurface,
BoundaryLayerSurface,
DiscreteSurface,
diff --git a/Geo/Makefile b/Geo/Makefile
index 6c85273c8d60bdb78424d68842e60506a1943e74..d2360106031956e6cb4fe546eb757b29f4766bf7 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.149 2007-07-09 13:54:36 geuzaine Exp $
+# $Id: Makefile,v 1.150 2007-07-26 16:28:27 anand Exp $
#
# Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
#
@@ -32,7 +32,7 @@ SRC = GEntity.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\
+ FProjectionFace.cpp\
GModel.cpp\
GModelIO_Geo.cpp\
GModelIO_Mesh.cpp\
@@ -70,213 +70,3 @@ depend:
rm -f Makefile.new
# DO NOT DELETE THIS LINE
-GEntity.o: GEntity.cpp GEntity.h Range.h SPoint3.h SBoundingBox3d.h \
- ../Common/GmshDefines.h MRep.h GEdge.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 GFace.h GEdgeLoop.h Pair.h \
- GRegion.h ../Common/VertexArray.h ../Common/Message.h ../Common/OS.h
-GVertex.o: GVertex.cpp GVertex.h GEntity.h Range.h SPoint3.h \
- SBoundingBox3d.h ../Common/GmshDefines.h MVertex.h GPoint.h SPoint2.h \
- GFace.h GEdgeLoop.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 Pair.h
-GEdge.o: GEdge.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
-GEdgeLoop.o: GEdgeLoop.cpp GEdgeLoop.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 ../Common/Message.h
-GFace.o: GFace.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
-GRegion.o: GRegion.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
-gmshVertex.o: gmshVertex.cpp 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 MEdge.h ../Common/Hash.h \
- MFace.h ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
- ExtrudeParams.h ../Common/SmoothData.h Pair.h gmshVertex.h Geo.h \
- gmshSurface.h ../DataStr/Tree.h ../DataStr/avl.h GeoInterpolation.h \
- ../Common/Message.h
-gmshEdge.o: gmshEdge.cpp 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 MEdge.h ../Common/Hash.h \
- MFace.h ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
- ExtrudeParams.h ../Common/SmoothData.h Pair.h gmshEdge.h Geo.h \
- gmshSurface.h ../DataStr/Tree.h ../DataStr/avl.h gmshVertex.h \
- GeoInterpolation.h ../Common/Message.h
-gmshFace.o: gmshFace.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 gmshVertex.h Geo.h gmshSurface.h ../DataStr/Tree.h \
- ../DataStr/avl.h gmshEdge.h gmshFace.h GeoInterpolation.h \
- ../Common/Message.h
-gmshRegion.o: gmshRegion.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 gmshFace.h Geo.h gmshSurface.h ../DataStr/Tree.h \
- ../DataStr/avl.h gmshVertex.h gmshRegion.h ../Common/Message.h
-gmshSurface.o: gmshSurface.cpp gmshSurface.h Pair.h Range.h SPoint2.h \
- SPoint3.h SVector3.h SBoundingBox3d.h ../Common/Message.h
-OCCVertex.o: OCCVertex.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 OCCVertex.h OCCIncludes.h OCCEdge.h OCCFace.h
-OCCEdge.o: OCCEdge.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 OCCEdge.h OCCVertex.h OCCIncludes.h \
- OCCFace.h
-OCCFace.o: OCCFace.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 OCCVertex.h OCCIncludes.h OCCEdge.h OCCFace.h \
- ../Common/Message.h
-OCCRegion.o: OCCRegion.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 OCCVertex.h OCCIncludes.h OCCEdge.h OCCFace.h OCCRegion.h \
- ../Common/Message.h
-FEdge.o: FEdge.cpp ../Common/Message.h 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 ../Common/Message.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 \
- MEdge.h ../Common/Hash.h MFace.h ../Numeric/Numeric.h \
- ../Common/Context.h ../DataStr/List.h ExtrudeParams.h \
- ../Common/SmoothData.h Pair.h
-GModel.o: GModel.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 gmshSurface.h ../Mesh/Field.h ../Post/Views.h \
- ../Post/ColorTable.h ../Common/VertexArray.h ../Post/AdaptiveViews.h \
- ../Common/GmshMatrix.h ../Geo/Geo.h ../Geo/gmshSurface.h \
- ../DataStr/Tree.h ../DataStr/avl.h ../Geo/SPoint2.h \
- ../Geo/ExtrudeParams.h ../Geo/GEdge.h ../Post/OctreePost.h \
- ../Common/Octree.h ../Common/OctreeInternals.h MRep.h \
- ../Common/Message.h ../Common/OS.h ../Mesh/BackgroundMesh.h
-GModelIO_Geo.o: GModelIO_Geo.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 Geo.h gmshSurface.h ../DataStr/Tree.h ../DataStr/avl.h \
- ../Parser/OpenFile.h ../DataStr/Tools.h ../DataStr/List.h \
- ../DataStr/Tree.h ../Common/Message.h gmshVertex.h gmshFace.h \
- gmshEdge.h gmshRegion.h ../Parser/Parser.h
-GModelIO_Mesh.o: GModelIO_Mesh.cpp ../Common/Message.h \
- ../Common/GmshDefines.h GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
- SBoundingBox3d.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 \
- gmshRegion.h Geo.h gmshSurface.h ../DataStr/Tree.h ../DataStr/avl.h \
- gmshFace.h gmshVertex.h gmshEdge.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 \
- ../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 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 ../Post/Views.h ../Post/ColorTable.h \
- ../Common/VertexArray.h ../Post/AdaptiveViews.h ../Common/GmshMatrix.h \
- FFace.h FEdge.h FVertex.h ../Mesh/meshGFace.h GModelIO_F.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 \
- 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 MNeighbour.h
-GModelIO_MED.o: GModelIO_MED.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
-ExtrudeParams.o: ExtrudeParams.cpp ../Common/Gmsh.h ../Common/Message.h \
- ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
- ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
- Geo.h ../Common/GmshDefines.h gmshSurface.h Pair.h Range.h SPoint2.h \
- SPoint3.h SVector3.h SBoundingBox3d.h ExtrudeParams.h \
- ../Common/SmoothData.h ../Numeric/Numeric.h
-Geo.o: Geo.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
- ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
- ../DataStr/List.h ../DataStr/Tree.h ../Numeric/Numeric.h Geo.h \
- ../Common/GmshDefines.h gmshSurface.h Pair.h Range.h SPoint2.h \
- SPoint3.h SVector3.h SBoundingBox3d.h ExtrudeParams.h \
- ../Common/SmoothData.h GModel.h GVertex.h GEntity.h MVertex.h GPoint.h \
- GEdge.h MElement.h MEdge.h ../Common/Hash.h MFace.h ../Common/Context.h \
- GFace.h GEdgeLoop.h GRegion.h GeoInterpolation.h
-GeoStringInterface.o: GeoStringInterface.cpp ../Common/Gmsh.h \
- ../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
- ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h \
- ../DataStr/Tree.h ../Numeric/Numeric.h Geo.h ../Common/GmshDefines.h \
- gmshSurface.h Pair.h Range.h SPoint2.h SPoint3.h SVector3.h \
- SBoundingBox3d.h ExtrudeParams.h ../Common/SmoothData.h \
- GeoStringInterface.h ../Parser/Parser.h ../Parser/OpenFile.h \
- ../Common/Context.h GModel.h GVertex.h GEntity.h MVertex.h GPoint.h \
- GEdge.h MElement.h MEdge.h ../Common/Hash.h MFace.h GFace.h GEdgeLoop.h \
- GRegion.h
-GeoInterpolation.o: GeoInterpolation.cpp ../Common/Gmsh.h \
- ../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
- ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h \
- ../DataStr/Tree.h Geo.h ../Common/GmshDefines.h gmshSurface.h Pair.h \
- Range.h SPoint2.h SPoint3.h SVector3.h SBoundingBox3d.h ExtrudeParams.h \
- ../Common/SmoothData.h ../Numeric/Numeric.h GeoInterpolation.h \
- GeoStringInterface.h
-findLinks.o: findLinks.cpp ../Common/Gmsh.h ../Common/Message.h \
- ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
- ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.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 ExtrudeParams.h ../Common/SmoothData.h GFace.h \
- GEdgeLoop.h Pair.h GRegion.h
-MVertex.o: MVertex.cpp MVertex.h SPoint3.h
-MElement.o: MElement.cpp MElement.h ../Common/GmshDefines.h MVertex.h \
- SPoint3.h MEdge.h SVector3.h ../Common/Hash.h MFace.h \
- ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h GEntity.h \
- Range.h SBoundingBox3d.h ../Common/Message.h
diff --git a/Geo/projectionFace.cpp b/Geo/projectionFace.cpp
deleted file mode 100644
index 5add1da1d1d3c65e1f62881a0118b37dd489e35e..0000000000000000000000000000000000000000
--- a/Geo/projectionFace.cpp
+++ /dev/null
@@ -1,180 +0,0 @@
-#include "projectionFace.h"
-
-SVector3 scalePoint(SVector3 p, SVector3 s)
-{
- double x = p[0]*s[0];
- double y = p[1]*s[1];
- double z = p[2]*s[2];
- SVector3 result(x,y,z);
- return result;
-}
-
-SVector3 rotatePoint(SVector3 p, SVector3 r)
-{
- double rot[16];
- double x = r[0] * Pi / 180.;
- double y = r[1] * Pi / 180.;
- double z = r[2] * Pi / 180.;
- double A = cos(x);
- double B = sin(x);
- double C = cos(y);
- double D = sin(y);
- double E = cos(z);
- double F = sin(z);
- double AD = A * D;
- double BD = B * D;
- rot[0] = C*E; rot[1] = BD*E+A*F; rot[2] =-AD*E+B*F; rot[3] = 0.;
- rot[4] =-C*F; rot[5] =-BD*F+A*E; rot[6] = AD*F+B*E; rot[7] = 0.;
- rot[8] = D; rot[9] =-B*C; rot[10] = A*C; rot[11] = 0.;
- rot[12] = 0.; rot[13] = 0.; rot[14] = 0.; rot[15] = 1.;
-
- x = p[0]*rot[0] + p[1]*rot[1] + p[2]*rot[2];
- y = p[0]*rot[4] + p[1]*rot[5] + p[2]*rot[6];
- z = p[0]*rot[8] + p[1]*rot[9] + p[2]*rot[10];
-
- SVector3 result(x,y,z);
- return result;
-}
-
-void projectionFace::rotate(SVector3 rot)
-{
- rotation = rot + rotation;
-}
-
-void projectionFace::translate(SVector3 trans)
-{
- translation = trans + translation;
-}
-
-void projectionFace::scale(SVector3 sc)
-{
- scaleFactor[0] *= sc[0];
- scaleFactor[1] *= sc[1];
- scaleFactor[2] *= sc[2];
-}
-
-void projectionFace::setScale(SVector3 sc)
-{
- scaleFactor[0] = sc[0];
- scaleFactor[1] = sc[1];
- scaleFactor[2] = sc[2];
-}
-
-void projectionFace::setRotation(SVector3 r)
-{
- rotation[0] = r[0];
- rotation[1] = r[1];
- rotation[2] = r[2];
-}
-
-void projectionFace::setTranslation(SVector3 t)
-{
- translation = t;
-}
-
-projectionFace::projectionFace(GModel *m,int num) : GFace(m,num)
-{
- for(int j = 0; j<3; j++){
- translation[j] = 0;
- rotation[j] = 0;
- }
-}
-
-projectionFace::~projectionFace()
-{
-}
-
-parabolicCylinder::parabolicCylinder(GModel *m, int num) : projectionFace(m,num)
-{
- focalPoint = 1;
- scaleFactor[0] = 1;
- scaleFactor[1] = 1;
- scaleFactor[2] = 1;
-}
-
-parabolicCylinder::~parabolicCylinder()
-{
-}
-
-GPoint parabolicCylinder::point(double par1, double par2) const
-{
- //first let's find the vertex
- //SPoint3 k = focalPoint - directrixPoint;
- //k.setPosition(k.x/2,k.y/2,k.z/2);
- //double f = sqrt(k.x^2 + k.y^2 + k.z^2);
- //SPoint3 vertex = directrixPoint + k;
- par1 -= .5; //I need these to make sure the 'vertex' occurs at the parameter point (.5,.5)
- par2 -= .5;
-
- SVector3 p(par1, (1/(4*focalPoint))*pow(par1,2.0), par2);
-
- //now we have to do the necessary transformations
- p = scalePoint(p,scaleFactor); //LOCAL SCALING
- p = rotatePoint(p,rotation);
- p = translation + p;
-
- GPoint gp(p[0],p[1],p[2]);
- return gp;
-}
-
-SPoint2 parabolicCylinder::parFromPoint(const SPoint3 &p) const
-{
- // ok...first we need to untranslate, unrotate and unscale it
- SVector3 trans(-translation[0], -translation[1], -translation[2]);
- SVector3 rot(-rotation[0],-rotation[1],-rotation[2]);
- SVector3 scalar(1/scaleFactor[0],1/scaleFactor[1],1/scaleFactor[2]); //LOCAL SCALING?
-
- SVector3 pos(p[0],p[1],p[2]);
- pos = trans + pos;
- pos = rotatePoint(pos,rot);
- pos = scalePoint(pos, scalar);
-
- SPoint2 q(pos[0],pos[1]);
- return q;
-}
-
-// du/dx, du/dy, du/dz; dv/dx, dv/dy, dv/dz
-// v - (1/4f)*u^2 = 0
-
-Pair<SVector3,SVector3> parabolicCylinder::firstDer(const SPoint2 ¶m) const
-{
- SVector3 du;
- SVector3 dv;
-
- du[0] = 1;
- du[1] = -(1/(2*focalPoint))*(param[0]-.5);
- du[2] = 0;
-
- dv[0] = 0;
- dv[1] = 0;
- dv[2] = 1;
-
- //ok, now we need to scale and rotate...hmm
- du = scalePoint(du, scaleFactor);
- dv = scalePoint(dv, scaleFactor);
- du = rotatePoint(du, rotation);
- dv = rotatePoint(dv, rotation);
- //yeah...that "might" work
-
- Pair<SVector3,SVector3> result(du,dv);
- return result;
-}
-
-// By convention, we're going to return the normal facing the interior
-// of the parabola
-SVector3 parabolicCylinder::normal(const SPoint2 ¶m) const
-{
- SVector3 n;
- Pair<SVector3,SVector3> result = firstDer(param);
-
- const SVector3 left = result.left();
- const SVector3 right = result.right();
- n = crossprod(left, right);
- n.normalize();
- return n;
-}
-
-Range<double> parabolicCylinder::parBounds(int i) const
-{
- return Range<double>(0, 1);
-}
diff --git a/Geo/projectionFace.h b/Geo/projectionFace.h
deleted file mode 100644
index 251bec9da07373f0c8aebae5688f7574e63a2460..0000000000000000000000000000000000000000
--- a/Geo/projectionFace.h
+++ /dev/null
@@ -1,62 +0,0 @@
-#ifndef _PROJECTION_FACE_H_
-#define _PROJECTION_FACE_H_
-
-#include "GFace.h"
-
-/* Specific to the derived classes, we're going to need functions
- which can be used to easily adjust the parameters of the functions.
- plus, we're going to have to implement all that stuff in GFace
- eventually. */
-SVector3 scalePoint(SVector3 p, SVector3 s);
-SVector3 rotatePoint(SVector3 p, SVector3 r);
-
-class projectionFace : public GFace {
- protected:
- SVector3 rotation; //this vector holds the euler angles at which the surface is rotated
- SVector3 translation; //this vector holds the location of the reference point in xyz space
- SVector3 scaleFactor; //this vector holds the scaling factors w.r.t x,y,z
- public:
- projectionFace(GModel *m, int num);
- ~projectionFace();
-
- void rotate(SVector3 rot); //rotates the surface by the euler angles rot
- void translate(SVector3 trans); //translates the surface by trans
- void scale(SVector3 sc); //scales the surface along the (local?) x,y,z axes by sc
- void setScale(SVector3 sc);
- void setRotation(SVector3 r);
- void setTranslation(SVector3 t);
-
- Range<double> parBounds(int i) const {throw;}
- virtual int paramDegeneracies(int dir, double *par) { return 0; }
- virtual GPoint point(double par1, double par2) const {throw;}
- virtual GPoint closestPoint(const SPoint3 & queryPoint) const {throw;}
- virtual int containsPoint(const SPoint3 &pt) const {throw;}
- virtual int containsParam(const SPoint2 &pt) const {throw;}
- virtual SVector3 normal(const SPoint2 ¶m) const {throw;}
- virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const {throw;}
- virtual GEntity::GeomType geomType() const { return GEntity::ProjectionSurface; }
- 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 UnknownModel; }
- void * getNativePtr() const {throw;}
- virtual bool surfPeriodic(int dim) const {throw;}
- virtual SPoint2 parFromPoint(const SPoint3 &) const {throw;}
-};
-
-class parabolicCylinder : public projectionFace {
- protected:
- double focalPoint; //the length from the vertex to the focal point
- public:
- parabolicCylinder(GModel *m, int num);
- ~parabolicCylinder();
- Range<double> parBounds(int i) const;
- GPoint point(double par1, double par2) const;
- SVector3 normal(const SPoint2 ¶m) const;
- Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
- SPoint2 parFromPoint(const SPoint3 &) const;
-};
-
-#endif
diff --git a/Graphics/Geom.cpp b/Graphics/Geom.cpp
index 9af8287a10317241a9aee58a0d02960db980bdfd..dfe90a6d644801a8d70cf229623c2349a594b4f0 100644
--- a/Graphics/Geom.cpp
+++ b/Graphics/Geom.cpp
@@ -1,4 +1,4 @@
-// $Id: Geom.cpp,v 1.130 2007-02-26 08:25:38 geuzaine Exp $
+// $Id: Geom.cpp,v 1.131 2007-07-26 16:28:27 anand Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -408,7 +408,7 @@ public :
if(f->geomType() == GEntity::Plane)
_drawPlaneGFace(f);
- else if(f->geomType() == GEntity::ProjectionSurface)
+ else if(f->geomType() == GEntity::ProjectionFace)
_drawProjectionGFace(f);
else
_drawNonPlaneGFace(f);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 87fb78f7514996ff9fe11107c6afdf3befa0732b..2d61c1af1a9e78ed323ad14339f71d72e4d9df13 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.81 2007-07-22 15:47:46 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.82 2007-07-26 16:28:27 anand Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -811,17 +811,17 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
// and compute characteristic lenghts using mesh edge spacing
BDS_GeomEntity CLASS_F (1,2);
-
+
it = edges.begin();
while(it != edges.end())
{
- if(!(*it)->is_mesh_degenerated()){
- if (!recover_medge ( m, *it))
- {
- Msg(GERROR,"Face not meshed");
- return false;
- }
- }
+ if(!(*it)->is_mesh_degenerated()){
+ if (!recover_medge ( m, *it))
+ {
+ Msg(GERROR,"Face not meshed");
+ return false;
+ }
+ }
++it;
}
// Msg(INFO,"Boundary Edges recovered for surface %d",gf->tag());
@@ -1530,7 +1530,7 @@ void meshGFace::operator() (GFace *gf)
{
if(gf->geomType() == GEntity::DiscreteSurface) return;
if(gf->geomType() == GEntity::BoundaryLayerSurface) return;
- if(gf->geomType() == GEntity::ProjectionSurface) return;
+ if(gf->geomType() == GEntity::ProjectionFace) return;
// destroy the mesh if it exists
deMeshGFace dem;
@@ -1602,7 +1602,7 @@ bool shouldRevert(MEdge &reference, std::vector<T*> &elements)
void orientMeshGFace::operator()(GFace *gf)
{
- if(gf->geomType() == GEntity::ProjectionSurface) return;
+ if(gf->geomType() == GEntity::ProjectionFace) return;
// surface orientions in OCC are not consistent with the orientation
// of the bounding edges, so just leave them unchanged:
diff --git a/contrib/FourierModel/BlendOperator.cpp b/contrib/FourierModel/BlendOperator.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..acaafd96013a5e06e53a579c9f8b92c538b0f1df
--- /dev/null
+++ b/contrib/FourierModel/BlendOperator.cpp
@@ -0,0 +1,451 @@
+#include "BlendOperator.h"
+
+bool BlendOperator::GetPointOnPatch_
+(int i, double u, double v, int j, double &x, double &y, double &z)
+{
+ double h = 1.e-8;
+ double tol = 1.e-12;
+
+ u = GetPatch(i)->RescaleU(u);
+ v = GetPatch(i)->RescaleV(v);
+
+ double P[3], d[3];
+ GetPatch(i)->GetProjectionSurface()->F(u,v,P[0],P[1],P[2]);
+ GetPatch(i)->GetProjectionSurface()->GetUnitNormal(u,v,d[0],d[1],d[2]);
+
+ u = GetPatch(i)->UnscaleU(u);
+ v = GetPatch(i)->UnscaleV(v);
+
+ double e1[3],e2[3];
+
+ if (std::abs(1.-std::abs(d[0])) < tol) {
+ e1[0] = 0.;
+ e1[1] = 1.;
+ e1[2] = 0.;
+ }
+ else if (std::abs(1.-std::abs(d[1])) < tol) {
+ e1[0] = 0.;
+ e1[1] = 0.;
+ e1[2] = 1.;
+ }
+ else if (std::abs(1.-std::abs(d[2])) < tol) {
+ e1[0] = 1.;
+ e1[1] = 0.;
+ e1[2] = 0.;
+ }
+ else {
+ e1[0] = - d[1] * d[2];
+ e1[1] = 2 * d[2] * d[0];
+ e1[2] = - d[0] * d[1];
+ }
+ double value = 0.;
+ for (int l=0;l<3;l++)
+ value += e1[l] * e1[l];
+ value = sqrt(value);
+ for (int l=0;l<3;l++)
+ e1[l] /= value;
+
+ e2[0] = d[1] * e1[2] - d[2] * e1[1];
+ e2[1] = d[2] * e1[0] - d[0] * e1[2];
+ e2[2] = d[0] * e1[1] - d[1] * e1[0];
+
+ bool converged = false;
+
+ double R[3];
+ GetPatch(i)->F(u,v,R[0],R[1],R[2]);
+
+ double Qu,Qv;
+ GetPatch(j)->GetProjectionSurface()->
+ OrthoProjectionOnSurface(R[0],R[1],R[2],Qu,Qv);
+
+ Qu = GetPatch(j)->UnscaleU(Qu);
+ Qv = GetPatch(j)->UnscaleV(Qv);
+
+ double Q[3];
+ GetPatch(j)->F(Qu,Qv,Q[0],Q[1],Q[2]);
+
+ int maxTotalIteration = 10;
+ int totalIterationCounter = 0;
+
+ while ((!converged) && (totalIterationCounter < maxTotalIteration)) {
+ totalIterationCounter++;
+ double QminusP[3];
+ for (int l=0;l<3;l++)
+ QminusP[l] = Q[l] - P[l];
+
+ double f[2];
+ f[0] = (QminusP[0] * e1[0] + QminusP[1] * e1[1] + QminusP[2] * e1[2]);
+ f[1] = (QminusP[0] * e2[0] + QminusP[1] * e2[1] + QminusP[2] * e2[2]);
+
+ //printf("(%g,%g) :: f = (%g,%g)\n",u,v,f[0],f[1]);
+
+ if ((std::abs(f[0]) < tol) && (std::abs(f[1]) < tol)) {
+ //bool IsGood = true;
+
+ bool IsGood = GetPatch(i)->GetProjectionSurface()->
+ OrthoProjectionOnSurface(Q[0],Q[1],Q[2],u,v);
+ if (IsGood) {
+ double RR[3];
+ //printf("(u,v) = (%g,%g); (Qu,Qv) = (%g,%g)\n",u,v,Qu,Qv);
+ u = GetPatch(i)->UnscaleU(u);
+ v = GetPatch(i)->UnscaleV(v);
+ GetPatch(i)->F(u,v,RR[0],RR[1],RR[2]);
+ double dist = 0.;
+ for (int l=0;l<3;l++)
+ dist += (RR[l] - R[l]) * (RR[l] - R[l]);
+ dist = sqrt(dist);
+ if (dist > 1.e-2)
+ break;
+ else {
+ //printf("R = (%g,%g,%g); RR = (%g,%g,%g)\n\n",
+ // R[0],R[1],R[2],RR[0],RR[1],RR[2]);
+ converged = true;
+ x = Q[0]; y = Q[1]; z = Q[2];
+ }
+ converged = true;
+ x = Q[0]; y = Q[1]; z = Q[2];
+ }
+ else
+ break;
+ }
+ else {
+ double QuPlus = Qu + h;
+ double QvPlus = Qv + h;
+ if (GetPatch(j)->_PI->periodic[0] == 1)
+ QuPlus -= std::floor(QuPlus);
+ if (GetPatch(j)->_PI->periodic[1] == 1)
+ QvPlus -= std::floor(QvPlus);
+
+ double QplusU[3], QplusV[3];
+ GetPatch(j)->F(QuPlus,Qv,QplusU[0],QplusU[1],QplusU[2]);
+ GetPatch(j)->F(Qu,QvPlus,QplusV[0],QplusV[1],QplusV[2]);
+
+ double QplusUminusP[3], QplusVminusP[3];
+ for (int l=0;l<3;l++) {
+ QplusUminusP[l] = QplusU[l] - P[l];
+ QplusVminusP[l] = QplusV[l] - P[l];
+ }
+
+ double fuPlus[2], fvPlus[2];
+ fuPlus[0] = (QplusUminusP[0] * e1[0] + QplusUminusP[1] * e1[1] +
+ QplusUminusP[2] * e1[2]);
+ fuPlus[1] = (QplusUminusP[0] * e2[0] + QplusUminusP[1] * e2[1] +
+ QplusUminusP[2] * e2[2]);
+ fvPlus[0] = (QplusVminusP[0] * e1[0] + QplusVminusP[1] * e1[1] +
+ QplusVminusP[2] * e1[2]);
+ fvPlus[1] = (QplusVminusP[0] * e2[0] + QplusVminusP[1] * e2[1] +
+ QplusVminusP[2] * e2[2]);
+
+ double D[2][2];
+ D[0][0] = (fuPlus[0] - f[0]) / h;
+ D[0][1] = (fvPlus[0] - f[0]) / h;
+ D[1][0] = (fuPlus[1] - f[1]) / h;
+ D[1][1] = (fvPlus[1] - f[1]) / h;
+
+ double det = D[0][0] * D[1][1] - D[0][1] * D[1][0];
+ double update[2];
+ update[0] = (D[1][1] * f[0] - D[0][1] * f[1]) / det;
+ update[1] = (D[0][0] * f[1] - D[1][0] * f[0]) / det;
+
+ //printf("update = (%g,%g)\n",update[0],update[1]);
+
+ Qu -= update[0];
+ Qv -= update[1];
+
+ //printf("Q1 = (%g,%g)\n",Qu,Qv);
+
+ if (GetPatch(j)->_PI->periodic[0] == 1)
+ Qu -= std::floor(Qu);
+ if (GetPatch(j)->_PI->periodic[1] == 1)
+ Qv -= std::floor(Qv);
+
+ //printf("Q2 = (%g,%g)\n",Qu,Qv);
+
+ GetPatch(j)->F(Qu,Qv,Q[0],Q[1],Q[2]);
+ }
+ }
+
+ return converged;
+}
+
+bool BlendOperator::E
+(int i, int j, double u, double v, double &x, double &y, double &z)
+{
+ bool found = false;
+
+ if (DoPatchesOverlap(i,j)) {
+ GetPatch(i)->F(u,v,x,y,z);
+ //if (1) {
+ if (IsPointInIntersectionBoundingBox(i,j,x,y,z)) {
+ if (GetPatch(i)->GetProjectionSurface()->GetTag() ==
+ GetPatch(j)->GetProjectionSurface()->GetTag()) {
+ found = GetPatch(j)->Inverse(x,y,z,u,v);
+ }
+ else {
+ int tag = GetProjectionSurfaceTagForOverlap(i,j);
+ if (GetPatch(i)->GetProjectionSurface()->GetTag() == tag) {
+ found = GetPointOnPatch_(i,u,v,j,x,y,z);
+ if (!GetPatch(j)->Inverse(x,y,z,u,v))
+ found = false;
+ }
+ else if (GetPatch(j)->GetProjectionSurface()->GetTag() == tag) {
+ found = GetPatch(j)->GetProjectionSurface()->
+ OrthoProjectionOnSurface(x,y,z,u,v);
+ if (found) {
+ u = GetPatch(j)->UnscaleU(u);
+ v = GetPatch(j)->UnscaleV(v);
+ double tol = 1.e-12;
+ if ((u > - tol)&&(u < 1. + tol)&&(v > - tol)&&(v < 1. + tol)) {
+ found = true;
+ GetPatch(j)->F(u,v,x,y,z);
+ }
+ else
+ found = false;
+ }
+ }
+ else {
+ found = E(i,tag,u,v,x,y,z);
+ if (found) {
+ found = GetPatch(tag)->Inverse(x,y,z,u,v);
+ if (found)
+ found = E(tag,j,u,v,x,y,z);
+ }
+ }
+ }
+ }
+ }
+ return found;
+}
+
+double BlendOperator::GetBlendingPou
+(int patchTag, double u, double v)
+{
+ double x,y,z;
+ double pouNumerator = GetPatch(patchTag)->GetPou(u,v);
+
+ double pouDenominator = 0.;
+ for (int otherPatchTag = 0;otherPatchTag < patch_.size();otherPatchTag++) {
+ if (E(patchTag,otherPatchTag,u,v,x,y,z)) {
+ double uu,vv;
+ GetPatch(otherPatchTag)->Inverse(x,y,z,uu,vv);
+ if (GetPatch(otherPatchTag)->Inverse(x,y,z,uu,vv)) {
+ pouDenominator += GetPatch(otherPatchTag)->GetPou(uu,vv);
+ }
+ }
+ }
+ if (pouDenominator)
+ return pouNumerator / pouDenominator;
+ else {
+ printf("%d(%g,%g) :: %g %g\n",patchTag,u,v,pouNumerator,pouDenominator);
+ return 1.;
+ }
+}
+
+Patch* BlendOperator::GetPatch
+(int tag) {
+ for (int i=0;i<patch_.size();i++) {
+ if (patch_[i]->GetTag() == tag)
+ return patch_[i];
+ }
+}
+
+ProjectionSurface* BlendOperator::GetProjectionSurface
+(int tag) {
+ for (int i=0;i<ps_.size();i++) {
+ if (ps_[i]->GetTag() == tag)
+ return ps_[i];
+ }
+}
+
+bool BlendOperator::DoPatchesOverlap
+(int patchTag1, int patchTag2) {
+ bool result;
+ if (overlapChart_[patchTag1][patchTag2]->doesIntersect)
+ result = true;
+ else
+ result = false;
+
+ return result;
+}
+
+int BlendOperator::GetProjectionSurfaceTagForOverlap
+(int patchTag1, int patchTag2) {
+ return overlapChart_[patchTag1][patchTag2]->psTag;
+}
+
+bool BlendOperator::IsPointInIntersectionBoundingBox
+(int patchTag1, int patchTag2, double x, double y, double z) {
+ bool result = true;
+ double tol = 1.e-2;
+ double xMin = overlapChart_[patchTag1][patchTag2]->xMin - tol;
+ double xMax = overlapChart_[patchTag1][patchTag2]->xMax + tol;
+ double yMin = overlapChart_[patchTag1][patchTag2]->yMin - tol;
+ double yMax = overlapChart_[patchTag1][patchTag2]->yMax + tol;
+ double zMin = overlapChart_[patchTag1][patchTag2]->zMin - tol;
+ double zMax = overlapChart_[patchTag1][patchTag2]->zMax + tol;
+
+ //printf("%d %d :: %g %g %g :: %g %g %g :: %g %g %g\n",
+ // patchTag1,patchTag2, xMin,x,xMax,yMin,y,yMax,zMin,z,zMax);
+
+ if ((x < xMin)||(x > xMax)||(y < yMin)||(y > yMax)||(z < zMin)||(z > zMax))
+ result = false;
+
+ return result;
+}
+
+int BlendOperator::GetNumPatches
+() {
+ return patch_.size();
+}
+
+bool BlendOperator::GetPointOnPatch
+(int patchTag, double d[3], double &x, double &y, double &z)
+{
+ double h = 1.e-6;
+ double tol = 1.e-12;
+
+ double P[3], R[3];
+ P[0] = R[0] = x;
+ P[1] = R[1] = y;
+ P[2] = R[2] = z;
+
+ double e1[3],e2[3];
+
+ if (std::abs(1.-std::abs(d[0])) < tol) {
+ e1[0] = 0.;
+ e1[1] = 1.;
+ e1[2] = 0.;
+ }
+ else if (std::abs(1.-std::abs(d[1])) < tol) {
+ e1[0] = 0.;
+ e1[1] = 0.;
+ e1[2] = 1.;
+ }
+ else if (std::abs(1.-std::abs(d[2])) < tol) {
+ e1[0] = 1.;
+ e1[1] = 0.;
+ e1[2] = 0.;
+ }
+ else {
+ e1[0] = - d[1] * d[2];
+ e1[1] = 2 * d[2] * d[0];
+ e1[2] = - d[0] * d[1];
+ }
+ double value = 0.;
+ for (int l=0;l<3;l++)
+ value += e1[l] * e1[l];
+ value = sqrt(value);
+ for (int l=0;l<3;l++)
+ e1[l] /= value;
+
+ e2[0] = d[1] * e1[2] - d[2] * e1[1];
+ e2[1] = d[2] * e1[0] - d[0] * e1[2];
+ e2[2] = d[0] * e1[1] - d[1] * e1[0];
+
+ bool converged = false;
+
+ double Qu,Qv;
+ GetPatch(patchTag)->GetProjectionSurface()->
+ OrthoProjectionOnSurface(P[0],P[1],P[2],Qu,Qv);
+
+
+ //printf("Q1 = (%g,%g)\n",Qu,Qv);
+
+ Qu = GetPatch(patchTag)->UnscaleU(Qu);
+ Qv = GetPatch(patchTag)->UnscaleV(Qv);
+
+ if (Qu < 0.)
+ Qu = 0.;
+ if (Qu > 1.)
+ Qu = 1.;
+
+ if (Qv < 0.)
+ Qv = 0.;
+ if (Qv > 1.)
+ Qv = 1.;
+
+ double Q[3];
+ GetPatch(patchTag)->F(Qu,Qv,Q[0],Q[1],Q[2]);
+
+ int maxTotalIteration = 10;
+ int totalIterationCounter = 0;
+
+ while ((!converged) && (totalIterationCounter < maxTotalIteration)) {
+ totalIterationCounter++;
+ double QminusP[3];
+ for (int l=0;l<3;l++)
+ QminusP[l] = Q[l] - P[l];
+
+ double f[2];
+ f[0] = (QminusP[0] * e1[0] + QminusP[1] * e1[1] + QminusP[2] * e1[2]);
+ f[1] = (QminusP[0] * e2[0] + QminusP[1] * e2[1] + QminusP[2] * e2[2]);
+
+ //printf("Q2 = (%g,%g) -> %g -> %g\n",Qu,Qv,f[0],f[1]);
+
+ if ((std::abs(f[0]) < tol) && (std::abs(f[1]) < tol)) {
+ //printf("here too\n");
+ converged = true;
+ for (int l=0;l<3;l++)
+ R[l] = Q[l];
+ }
+ else {
+ double QuPlus = Qu + h;
+ double QvPlus = Qv + h;
+ if (GetPatch(patchTag)->_PI->periodic[0] == 1)
+ QuPlus -= std::floor(QuPlus);
+ if (GetPatch(patchTag)->_PI->periodic[1] == 1)
+ QvPlus -= std::floor(QvPlus);
+
+ double QplusU[3], QplusV[3];
+ GetPatch(patchTag)->F(QuPlus,Qv,QplusU[0],QplusU[1],QplusU[2]);
+ GetPatch(patchTag)->F(Qu,QvPlus,QplusV[0],QplusV[1],QplusV[2]);
+
+ double QplusUminusP[3], QplusVminusP[3];
+ for (int l=0;l<3;l++) {
+ QplusUminusP[l] = QplusU[l] - P[l];
+ QplusVminusP[l] = QplusV[l] - P[l];
+ }
+
+ double fuPlus[2], fvPlus[2];
+ fuPlus[0] = (QplusUminusP[0] * e1[0] + QplusUminusP[1] * e1[1] +
+ QplusUminusP[2] * e1[2]);
+ fuPlus[1] = (QplusUminusP[0] * e2[0] + QplusUminusP[1] * e2[1] +
+ QplusUminusP[2] * e2[2]);
+ fvPlus[0] = (QplusVminusP[0] * e1[0] + QplusVminusP[1] * e1[1] +
+ QplusVminusP[2] * e1[2]);
+ fvPlus[1] = (QplusVminusP[0] * e2[0] + QplusVminusP[1] * e2[1] +
+ QplusVminusP[2] * e2[2]);
+
+ double D[2][2];
+ D[0][0] = (fuPlus[0] - f[0]) / h;
+ D[0][1] = (fvPlus[0] - f[0]) / h;
+ D[1][0] = (fuPlus[1] - f[1]) / h;
+ D[1][1] = (fvPlus[1] - f[1]) / h;
+
+ double det = D[0][0] * D[1][1] - D[0][1] * D[1][0];
+ double update[2];
+ update[0] = (D[1][1] * f[0] - D[0][1] * f[1]) / det;
+ update[1] = (D[0][0] * f[1] - D[1][0] * f[0]) / det;
+
+ //printf("update = (%g,%g)\n",update[0],update[1]);
+
+ Qu -= update[0];
+ Qv -= update[1];
+
+ //printf("Q1 = (%g,%g)\n",Qu,Qv);
+
+ //if (GetPatch(patchTag)->_PI->periodic[0] == 1)
+ //Qu -= std::floor(Qu);
+ //if (GetPatch(patchTag)->_PI->periodic[1] == 1)
+ //Qv -= std::floor(Qv);
+
+ //printf("Q2 = (%g,%g)\n",Qu,Qv);
+
+ GetPatch(patchTag)->F(Qu,Qv,Q[0],Q[1],Q[2]);
+ }
+ }
+ x = R[0];
+ y = R[1];
+ z = R[2];
+ return converged;
+}
diff --git a/contrib/FourierModel/BlendOperator.h b/contrib/FourierModel/BlendOperator.h
new file mode 100755
index 0000000000000000000000000000000000000000..db8b1f593a254c4687102b0a26b942ab5df13401
--- /dev/null
+++ b/contrib/FourierModel/BlendOperator.h
@@ -0,0 +1,56 @@
+#ifndef _BLEND_OPERATOR_H_
+#define _BLEND_OPERATOR_H_
+
+#include <vector>
+#include "FM_Info.h"
+#include "Patch.h"
+#include "ProjectionSurface.h"
+
+class BlendOperator {
+ private:
+ std::vector<Patch*> patch_;
+ std::vector<ProjectionSurface*> ps_;
+ std::vector<std::vector<OverlapInfo*> > overlapChart_;
+
+ bool GetPointOnPatch_
+ (int i, double u, double v, int j, double &x, double &y, double &z);
+
+ protected:
+
+ public:
+ BlendOperator
+ (std::vector<Patch*> patch,std::vector<ProjectionSurface*> ps,
+ std::vector<std::vector<OverlapInfo*> > overlapChart) :
+ patch_(patch), ps_(ps), overlapChart_(overlapChart) {}
+
+ ~BlendOperator() {}
+
+ int GetNumPatches
+ ();
+
+ Patch* GetPatch
+ (int tag);
+
+ ProjectionSurface* GetProjectionSurface
+ (int tag);
+
+ bool DoPatchesOverlap
+ (int patchTag1, int patchTag2);
+
+ int GetProjectionSurfaceTagForOverlap
+ (int patchTag1, int patchTag2);
+
+ bool IsPointInIntersectionBoundingBox
+ (int patchTag1, int patchTag2, double x, double y, double z);
+
+ bool E
+ (int i, int j, double u, double v, double &x, double &y, double &z);
+
+ double GetBlendingPou
+ (int patchTag, double u, double v);
+
+ bool GetPointOnPatch
+ (int patchTag, double d[3], double &x, double &y, double &z);
+};
+
+#endif
diff --git a/contrib/FourierModel/BlendedPatch.cpp b/contrib/FourierModel/BlendedPatch.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..ecfeb0e323c64d4673239b3d04fa1e442c1309aa
--- /dev/null
+++ b/contrib/FourierModel/BlendedPatch.cpp
@@ -0,0 +1,878 @@
+#include "BlendedPatch.h"
+
+BlendedPatch::BlendedPatch
+(Patch* patch, BlendOperator* blendOp)
+ : _coeffData(0),_coeffDerivU(0),_coeffDerivV(0),
+ _coeffDerivUU(0),_coeffDerivVV(0),_coeffDerivUV(0)
+{
+ patch_ = patch;
+ blendOp_ = blendOp;
+
+ patchTag_ = patch_->GetTag();
+ nPatches_ = blendOp_->GetNumPatches();
+
+ _derivative = patch_->GetDerivativeBitField();
+
+ _uM = patch_->GetUModes();
+ _vM = patch_->GetVModes();
+
+ _periodU = patch_->IsUPeriodic() ? 1. : 2.;
+ _periodV = patch_->IsVPeriodic() ? 1. : 2.;
+
+ IsUPeriodic = false;
+ IsVPeriodic = false;
+
+ _uM = 32;
+ _vM = 32;
+
+ //printf("%d %d\n",_uM,_vM);
+
+ if (IsUPeriodic) {
+ if ((_uM % 2) == 0) {
+ _uMLower = - _uM/2;
+ _uMUpper = _uM/2 - 1;
+ }
+ else {
+ _uMLower = - (_uM - 1)/2;
+ _uMUpper = (_uM - 1)/2;
+ }
+ }
+ else {
+ _uMLower = 0;
+ _uMUpper = _uM;
+ }
+ if (IsVPeriodic) {
+ if ((_vM % 2) == 0) {
+ _vMLower = - _vM/2;
+ _vMUpper = _vM/2 - 1;
+ }
+ else {
+ _vMLower = - (_vM - 1)/2;
+ _vMUpper = (_vM - 1)/2;
+ }
+ }
+ else {
+ _vMLower = 0;
+ _vMUpper = _vM;
+ }
+
+ if (_derivative)
+ _ProcessSeriesCoeff();
+
+ /*
+ for (int j = 0; j < _uM; j++) {
+ for (int k = 0; k < _vM; k++)
+ printf("%g ",_coeffData[j][k].real());
+ printf("\n\n");
+ }
+ printf("\n");
+ */
+
+ // 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>*[_uM];
+ _coeffDerivV = new std::complex<double>*[_uM];
+ for(int j = 0; j < _uM; j++){
+ _coeffDerivU[j] = new std::complex<double>[_vM];
+ _coeffDerivV[j] = new std::complex<double>[_vM];
+ for(int k = 0; k < _vM; k++){
+ _coeffDerivU[j][k] = 0.;
+ _coeffDerivV[j][k] = 0.;
+ }
+ }
+ }
+
+ if(_derivative & 2){
+ _coeffDerivUU = new std::complex<double>*[_uM];
+ _coeffDerivVV = new std::complex<double>*[_uM];
+ _coeffDerivUV = new std::complex<double>*[_uM];
+ for(int j = 0; j < _uM; j++){
+ _coeffDerivUU[j] = new std::complex<double>[_vM];
+ _coeffDerivVV[j] = new std::complex<double>[_vM];
+ _coeffDerivUV[j] = new std::complex<double>[_vM];
+ for(int k = 0; k < _vM; 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 = _uM - 1; j >= 0; j--){
+ for(int k = _vM - 1; k >= 0; k--){
+ if(_derivative & 1){
+ if (IsUPeriodic) {
+ int J = j+_uMLower;
+ _coeffDerivU[j][k] = (2 * M_PI * J * I / _periodU) *
+ _coeffData[j][k];
+ }
+ else {
+ if (j == _uM - 1)
+ _coeffDerivU[j][k] = 0.;
+ else if (j == _uM - 2)
+ _coeffDerivU[j][k] = 2. * (double)(j + 1) * _coeffData[j + 1][k];
+ else
+ _coeffDerivU[j][k] = _coeffDerivU[j + 2][k] +
+ 2. * (double)(j + 1) * _coeffData[j + 1][k];
+ //if (j != 0)
+ //_coeffDerivU[j][k] *= 2.;
+ }
+ if (IsVPeriodic) {
+ int K = k+_vMLower;
+ _coeffDerivV[j][k] = (2 * M_PI * K * I / _periodV) *
+ _coeffData[j][k];
+ }
+ else {
+ if (k == _vM - 1)
+ _coeffDerivV[j][k] = 0.;
+ else if (k == _vM - 2)
+ _coeffDerivV[j][k] = 2. * (double)(k + 1) * _coeffData[j][k + 1];
+ else
+ _coeffDerivV[j][k] = _coeffDerivV[j][k + 2] +
+ 2. * (double)(k + 1) * _coeffData[j][k + 1];
+ //if (k != 0)
+ //_coeffDerivV[j][k] *= 2.;
+ }
+ }
+ }
+ }
+ for (int j = 0; j < _uM; j++) {
+ for (int k = 0; k < _vM; k++) {
+ if (_derivative & 1) {
+ if (!IsUPeriodic) {
+ if (j != 0) {
+ _coeffDerivU[j][k] *= 2.;
+ }
+ }
+ if (!IsVPeriodic) {
+ if (k != 0) {
+ _coeffDerivV[j][k] *= 2.;
+ }
+ }
+ }
+ }
+ }
+ for(int j = _uM - 1; j >= 0; j--) {
+ for(int k = _vM - 1; k >= 0; k--) {
+ if(_derivative & 2) {
+ if (IsUPeriodic) {
+ int J = j+_uMLower;
+ _coeffDerivUU[j][k] = (2 * M_PI * J * I / _periodU) *
+ _coeffDerivU[j][k];
+ }
+ else {
+ if (j == _uM - 1)
+ _coeffDerivUU[j][k] = 0.;
+ else if (j == _uM - 2)
+ _coeffDerivUU[j][k] = 2. * (double)(j + 1) *
+ _coeffDerivU[j + 1][k];
+ else
+ _coeffDerivUU[j][k] = _coeffDerivUU[j + 2][k] +
+ 2. * (double)(j + 1) * _coeffDerivU[j + 1][k];
+ //if (j != 0)
+ //_coeffDerivUU[j][k] *= 2.;
+ }
+ if (IsVPeriodic) {
+ int K = k+_vMLower;
+ _coeffDerivVV[j][k] = (2 * M_PI * K * I / _periodV) *
+ _coeffDerivV[j][k];
+ _coeffDerivUV[j][k] = (2 * M_PI * K * I / _periodV) *
+ _coeffDerivU[j][k];
+ }
+ else {
+ if (k == _vM - 1) {
+ _coeffDerivVV[j][k] = 0.;
+ _coeffDerivUV[j][k] = 0.;
+ }
+ else if (k == _vM - 2) {
+ _coeffDerivVV[j][k] = 2. * (double)(k + 1) *
+ _coeffDerivV[j][k + 1];
+ _coeffDerivUV[j][k] = 2. * (double)(k + 1) *
+ _coeffDerivU[j][k + 1];
+ }
+ else {
+ _coeffDerivVV[j][k] = _coeffDerivVV[j][k + 2] +
+ 2. * (double)(k + 1) * _coeffDerivV[j][k + 1];
+ _coeffDerivUV[j][k] = _coeffDerivUV[j][k + 2] +
+ 2. * (double)(k + 1) * _coeffDerivU[j][k + 1];
+ //if (k != 0) {
+ //_coeffDerivVV[j][k] *= 2.;
+ //_coeffDerivUV[j][k] *= 2.;
+ //}
+ }
+ }
+ }
+ }
+ }
+ for (int j = 0; j < _uM; j++) {
+ for (int k = 0; k < _vM; k++) {
+ if (_derivative & 2) {
+ if (!IsUPeriodic && IsVPeriodic) {
+ if (j != 0) {
+ _coeffDerivUU[j][k] *= 2.;
+ }
+ }
+ if (IsUPeriodic && !IsVPeriodic) {
+ if (k != 0) {
+ _coeffDerivVV[j][k] *= 2.;
+ _coeffDerivUV[j][k] *= 2.;
+ }
+ }
+ }
+ }
+ }
+ }
+ // Initialize interpolation variables
+ _tmpCoeff = std::vector< std::complex<double> >(_vM);
+ _tmpInterp = std::vector< std::complex<double> >(_uM);
+}
+
+BlendedPatch::~BlendedPatch()
+{
+ for(int j = 0; j < _uM; j++)
+ delete [] _coeffData[j];
+ delete [] _coeffData;
+
+ if(_coeffDerivU){
+ for(int j = 0; j < _uM; j++)
+ delete [] _coeffDerivU[j];
+ delete [] _coeffDerivU;
+ }
+ if(_coeffDerivV){
+ for(int j = 0; j < _uM; j++)
+ delete [] _coeffDerivV[j];
+ delete [] _coeffDerivV;
+ }
+ if(_coeffDerivUU){
+ for(int j = 0; j < _uM; j++)
+ delete [] _coeffDerivUU[j];
+ delete [] _coeffDerivUU;
+ }
+ if(_coeffDerivVV){
+ for(int j = 0; j < _uM; j++)
+ delete [] _coeffDerivVV[j];
+ delete [] _coeffDerivVV;
+ }
+ if(_coeffDerivUV){
+ for(int j = 0; j < _uM; j++)
+ delete [] _coeffDerivUV[j];
+ delete [] _coeffDerivUV;
+ }
+}
+
+int BlendedPatch::_forwardSize = 0;
+int BlendedPatch::_backwardSize = 0;
+fftw_plan BlendedPatch::_forwardPlan;
+fftw_plan BlendedPatch::_backwardPlan;
+fftw_complex *BlendedPatch::_forwardData = 0;
+fftw_complex *BlendedPatch::_backwardData = 0;
+
+void BlendedPatch::_SetForwardPlan(int n)
+{
+ if(n != _forwardSize){
+ if(_forwardSize){
+ fftw_destroy_plan(_forwardPlan);
+ fftw_free(_forwardData);
+ }
+ _forwardSize = n;
+ _forwardData =
+ (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * _forwardSize);
+ _forwardPlan =
+ fftw_plan_dft_1d(_forwardSize, _forwardData, _forwardData,
+ FFTW_FORWARD, FFTW_ESTIMATE);
+ }
+}
+
+void BlendedPatch::_SetBackwardPlan(int n)
+{
+ if(n != _backwardSize){
+ if(_backwardSize){
+ fftw_destroy_plan(_backwardPlan);
+ fftw_free(_backwardData);
+ }
+ _backwardSize = n;
+ _backwardData =
+ (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * _backwardSize);
+ _backwardPlan =
+ fftw_plan_dft_1d(_backwardSize, _backwardData, _backwardData,
+ FFTW_BACKWARD, FFTW_ESTIMATE);
+ }
+}
+
+void BlendedPatch::_ForwardFft(int n, std::complex<double> *fftData)
+{
+ // Initialize fftw plan and array (ignoring the last element of
+ // fftData, which should just be the periodic extension)
+ _SetForwardPlan(n - 1);
+ for(int i = 0; i < n - 1; i++){
+ _forwardData[i][0] = fftData[i].real();
+ _forwardData[i][1] = fftData[i].imag();
+ }
+
+ // Perform forward FFT
+ fftw_execute(_forwardPlan);
+
+ // Copy data back into fftData and scale by 1/(n - 1)
+ double s = 1. / (double)(n - 1);
+ for(int i = 0; i < n - 1; i++)
+ fftData[i] =
+ s * std::complex<double>(_forwardData[i][0], _forwardData[i][1]);
+}
+
+void BlendedPatch::_BackwardFft(int n, std::complex<double> *fftData)
+{
+ // Initialize fftw plan and array (ignoring last element of fftData)
+ _SetBackwardPlan(n - 1);
+ for(int i = 0; i < n - 1; i++){
+ _backwardData[i][0] = fftData[i].real();
+ _backwardData[i][1] = fftData[i].imag();
+ }
+
+ // Perform backward FFT
+ fftw_execute(_backwardPlan);
+
+ // Copy data back into fftData
+ for(int i = 0; i < n - 1; i++)
+ fftData[i] =
+ std::complex<double>(_backwardData[i][0], _backwardData[i][1]);
+
+ // Fill in last element with copy of first element
+ fftData[n - 1] = fftData[0];
+}
+
+void BlendedPatch::_ProcessSeriesCoeff()
+{
+ _coeffData = new std::complex<double>*[_uM];
+ for(int j = 0; j < _uM; j++) {
+ _coeffData[j] = new std::complex<double>[_vM];
+ for(int k = 0; k < _vM; k++)
+ _coeffData[j][k] = 0.;
+ }
+ if (IsUPeriodic) {
+ std::vector<double> u(_uM), v(_vM + 1);
+ for (int j = 0; j < _uM; j++)
+ u[j] = (double)j / (double)(_uM-1);
+ for (int j = 0; j < _vM + 1; j++)
+ v[j] = (double)j / (double)_vM;
+
+ std::complex<double> **dataU = new std::complex<double> *[_vM];
+ for (int k = 0; k < _vM; k++)
+ dataU[k] = new std::complex<double> [_uM];
+ std::complex<double> *dataV = new std::complex<double>[2*_vM + 1];
+ for (int j = 0; j < _uM - 1; j++) {
+ for (int k = 0; k < _vM + 1; k++) {
+
+ double X,Y,Z,pX,pY,pZ,nX,nY,nZ;
+ F(u[j],0.5 * cos(M_PI * v[k]) + 0.5,X,Y,Z);
+ patch_->GetProjectionSurface()->
+ F(patch_->RescaleU(u[j]),
+ patch_->RescaleV(0.5 * cos(M_PI * v[k]) + 0.5),pX,pY,pZ);
+ patch_->GetProjectionSurface()->GetUnitNormal
+ (patch_->RescaleU(u[j]),
+ patch_->RescaleV(0.5 * cos(M_PI * v[k]) + 0.5),nX,nY,nZ);
+ dataV[k] = (X-pX)*nX + (Y-pY)*nY + (Z-pZ)*nZ;
+
+ //dataV[k] = (0.5 * cos(M_PI * u[j]) + 0.5) *
+ //(0.5 * cos(M_PI * v[k]) + 0.5);
+ }
+
+ for (int k = 1; k < _vM+1; k++)
+ dataV[_vM + k] = dataV[_vM -k];
+ _BackwardFft(2*_vM + 1, dataV);
+ dataU[0][j] = 0.5 * dataV[0] / (double)_vM;
+ for (int k=1; k<_vM-1; k++)
+ dataU[k][j] = dataV[k] / (double)_vM;
+ dataU[_vM-1][j] = 0.5 * dataV[_vM-1] / (double)_vM;
+ }
+
+ for (int k = 0; k < _vM; k++) {
+ dataU[k][_uM - 1] = dataU[k][0];
+ _ForwardFft(_uM, dataU[k]);
+ }
+
+ for (int j = _uMLower; j <= _uMUpper; j++) {
+ for (int k = 0; k < _vM; k++)
+ if ((j == _uMLower) || (j == _uMUpper))
+ _coeffData[_uMUpper + j][k] = dataU[k][_uMUpper] / 2.;
+ else if ((j >= 0) && (j < _uMUpper))
+ _coeffData[_uMUpper + j][k] = dataU[k][j];
+ else
+ _coeffData[_uMUpper + j][k] = dataU[k][_uM + j -1];
+ }
+ for (int k = 0; k < _vM; k++)
+ delete [] dataU[k];
+ delete [] dataU;
+ delete [] dataV;
+ }
+ else if (IsVPeriodic) {
+ std::vector<double> u(_uM + 1), v(_vM);
+ for (int j = 0; j < _uM + 1; j++)
+ u[j] = (double)j / (double)_uM;
+ for (int j = 0; j < _vM; j++)
+ v[j] = (double)j / (double)(_vM-1);
+
+ std::complex<double> **dataV = new std::complex<double> *[_uM];
+ for (int j = 0; j < _uM; j++)
+ dataV[j] = new std::complex<double> [_vM];
+ std::complex<double> *dataU = new std::complex<double>[2*_uM + 1];
+ for (int k = 0; k < _vM - 1; k++) {
+ for (int j = 0; j < _uM + 1; j++) {
+ double X,Y,Z,pX,pY,pZ,nX,nY,nZ;
+ F(0.5 * cos(M_PI * u[j]) + 0.5,v[k],X,Y,Z);
+ patch_->GetProjectionSurface()->
+ F(patch_->RescaleU(0.5 * cos(M_PI * u[j]) + 0.5),
+ patch_->RescaleV(v[k]),pX,pY,pZ);
+ patch_->GetProjectionSurface()->GetUnitNormal
+ (patch_->RescaleU(u[j]),
+ patch_->RescaleV(0.5 * cos(M_PI * v[k]) + 0.5),nX,nY,nZ);
+ dataU[j] = (X-pX)*nX + (Y-pY)*nY + (Z-pZ)*nZ;
+ }
+ for (int j = 1; j < _uM+1; j++)
+ dataU[_uM + j] = dataU[_uM -j];
+ _BackwardFft(2*_uM + 1, dataU);
+ dataV[0][k] = 0.5 * dataU[0] / (double)_uM;
+ for (int j=1; j<_uM-1; j++)
+ dataV[j][k] = dataU[j] / (double)_uM;
+ dataV[_uM-1][k] = 0.5 * dataU[_uM-1] / (double)_uM;
+ }
+ for (int j = 0; j < _uM; j++) {
+ dataV[j][_uM - 1] = dataV[j][0];
+ _ForwardFft(_vM, dataV[j]);
+ }
+ for (int k = _vMLower; k <= _vMUpper; k++) {
+ //for (int j = 0; j < _uM - 1; j++) {
+ for (int j = 0; j < _uM; j++)
+ if ((k == _vMLower) || (k == _vMUpper))
+ _coeffData[j][_vMUpper + k] = dataV[j][_vMUpper] / 2.;
+ else if ((j >= 0) && (j < _vMUpper))
+ _coeffData[j][_uMUpper + k] = dataV[j][k];
+ else
+ _coeffData[j][_uMUpper + k] = dataV[j][_uM + k -1];
+ }
+ for (int j = 0; j < _uM; j++)
+ delete [] dataV[j];
+ delete [] dataV;
+ delete [] dataU;
+ }
+ else {
+ std::vector<double> u(_uM + 1), v(_vM + 1);
+ for (int j = 0; j < _uM + 1; j++)
+ u[j] = (double)j / (double)_uM;
+ for (int j = 0; j < _vM + 1; j++)
+ v[j] = (double)j / (double)_vM;
+
+ std::complex<double> **dataU = new std::complex<double> *[_vM];
+ for (int k = 0; k < _vM; k++)
+ dataU[k] = new std::complex<double> [2*_uM + 1];
+
+ std::complex<double> *dataV = new std::complex<double>[2*_vM + 1];
+ for (int j = 0; j < _uM + 1; j++) {
+ for (int k = 0; k < _vM + 1; k++) {
+ /*
+ double X,Y,Z,pX,pY,pZ,nX,nY,nZ;
+ _F(0.5 * cos(M_PI * u[j]) + 0.5,0.5 * cos(M_PI * v[k]) + 0.5,X,Y,Z);
+ patch_->GetProjectionSurface()->
+ F(patch_->RescaleU(0.5 * cos(M_PI * u[j]) + 0.5),
+ patch_->RescaleV(0.5 * cos(M_PI * v[k]) + 0.5),pX,pY,pZ);
+ patch_->GetProjectionSurface()->GetUnitNormal
+ (patch_->RescaleU(0.5 * cos(M_PI * u[j]) + 0.5),
+ patch_->RescaleV(0.5 * cos(M_PI * v[k]) + 0.5),nX,nY,nZ);
+ dataV[k] = (X-pX)*nX + (Y-pY)*nY + (Z-pZ)*nZ;
+ */
+ //dataV[k] = (0.5 * cos(M_PI * u[j]) + 0.5) *
+ //(0.5 * cos(M_PI * v[k]) + 0.5);
+ dataV[k] = (0.5 * cos(M_PI * v[k]) + 0.5);
+ }
+ for (int k = 1; k < _vM+1; k++)
+ dataV[_vM + k] = dataV[_vM - k];
+ _BackwardFft(2*_vM + 1, dataV);
+ dataU[0][j] = 0.5 * dataV[0] / (double)_vM;
+ for (int k=1; k<_vM-1; k++)
+ dataU[k][j] = dataV[k] / (double)_vM;
+ dataU[_vM-1][j] = 0.5 * dataV[_vM-1] / (double)_vM;
+ }
+
+ /*
+ for (int k = 0; k < _vM; k++)
+ for (int j = 0; j < _uM + 1; j++)
+ printf("%d %d -> %g %g\n",k,j,dataU[k][j].real(),dataU[k][j].imag());
+ */
+
+ for (int k = 0; k < _vM; k++) {
+ for (int j = 1; j < _uM+1; j++) {
+ dataU[k][_uM + j] = dataU[k][_uM - k];
+ }
+ _BackwardFft(2*_uM + 1, dataU[k]);
+ }
+
+ for (int j = 0; j < _uM; j++) {
+ for (int k = 0; k < _vM; k++) {
+ if ((j == 0) || (j == _uM - 1))
+ _coeffData[j][k] = 0.5 * dataU[k][j]/ (double)_uM;
+ else
+ _coeffData[j][k] = dataU[k][j]/ (double)_uM;
+ }
+ }
+ for (int k = 0; k < _vM; k++)
+ delete [] dataU[k];
+ delete [] dataU;
+ delete [] dataV;
+ }
+}
+
+std::complex<double> BlendedPatch::
+_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> BlendedPatch::
+ _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.);
+ }
+ std::vector<double> uT(_uM,0.);
+ std::vector<double> vT(_vM,0.);
+ if (!IsUPeriodic) {
+ for (int j = 0; j < _uM; j++)
+ if (j == 0)
+ uT[j] = 1.;
+ else if (j == 1)
+ uT[j] = 2. * u - 1;
+ else
+ uT[j] = 2. * uT[1] * uT[j-1] - uT[j-2];
+ }
+ if (!IsVPeriodic) {
+ for (int k = 0; k < _vM; k++)
+ if (k == 0)
+ vT[k] = 1.;
+ else if (k == 1)
+ vT[k] = 2. * v - 1.;
+ else
+ vT[k] = 2. * vT[1] * vT[k-1] - vT[k-2];
+ }
+ // Interpolate to find value at (u,v)
+ for(int j = 0; j < _uM; j++){
+ _tmpInterp[j] = 0.;
+ for(int k = 0; k < _vM; 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;
+ }
+ if (IsVPeriodic) {
+ 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 * _vMLower * v / _periodV),
+ sin(2 * M_PI * _vMLower * v / _periodV));
+ }
+ else {
+ //printf("i was here 0\n");
+ for(int k = 0; k < _vM; k++)
+ _tmpInterp[j] += _tmpCoeff[k] * vT[k];
+ }
+ }
+ if (IsUPeriodic) {
+ 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 * _uMLower * u / _periodU),
+ sin(2 * M_PI * _uMLower * u / _periodU));
+ }
+ else {
+ std::complex<double> tmp = 0.;
+ for(int j = 0; j < _uM; j++)
+ tmp += _tmpInterp[j] * uT[j];
+ return tmp;
+ }
+}
+
+void BlendedPatch::F
+(double u, double v, double &x, double &y, double &z)
+{
+ x = y = z = 0.;
+
+ for (int otherPatchTag_ = 0;otherPatchTag_ < nPatches_; otherPatchTag_++) {
+ double xx, yy, zz;
+ if (blendOp_->E(patchTag_,otherPatchTag_,u,v,xx,yy,zz)) {
+ double uu,vv;
+ blendOp_->GetPatch(otherPatchTag_)->Inverse(xx,yy,zz,uu,vv);
+ double lambda = blendOp_->GetBlendingPou(otherPatchTag_,uu,vv);
+ x += lambda * xx;
+ y += lambda * yy;
+ z += lambda * zz;
+ }
+ }
+}
+
+bool BlendedPatch::Inverse
+(double x, double y, double z, double &u, double &v)
+{
+ bool found = false;
+ int projectionSurfaceTag_;
+ for (int otherPatchTag_ = 0;otherPatchTag_ < nPatches_; otherPatchTag_++) {
+ if (otherPatchTag_ != patchTag_) {
+ if (blendOp_->DoPatchesOverlap(patchTag_,otherPatchTag_)) {
+ if (blendOp_->
+ IsPointInIntersectionBoundingBox(patchTag_,otherPatchTag_,x,y,z)) {
+ projectionSurfaceTag_ = blendOp_->
+ GetProjectionSurfaceTagForOverlap(patchTag_,otherPatchTag_);
+ found = true;
+ break;
+ }
+ }
+ }
+ }
+ if (found) {
+ found = blendOp_->GetProjectionSurface(projectionSurfaceTag_)->
+ OrthoProjectionOnSurface(x,y,z,u,v);
+ double xx,yy,zz;
+ blendOp_->GetProjectionSurface(projectionSurfaceTag_)->F(u,v,xx,yy,zz);
+ double d[3];
+ d[0] = x - xx;
+ d[1] = y - yy;
+ d[2] = z - zz;
+ double abs_d = 0.;
+ for (int l=0;l<3;l++)
+ abs_d += d[l] * d[l];
+ abs_d = std::sqrt(abs_d);
+ for (int l=0;l<3;l++)
+ d[l] /= abs_d;
+ found = blendOp_->GetPointOnPatch(patchTag_,d,x,y,z);
+ }
+ found = blendOp_->GetPatch(patchTag_)->Inverse(x,y,z,u,v);
+}
+
+/*
+void BlendedPatch::F
+(double u, double v, double &x, double &y, double &z)
+{
+ if (_derivative) {
+ double px, py, pz, nx, ny, nz, d;
+
+ u = patch_->RescaleU(u);
+ v = patch_->RescaleV(v);
+
+ patch_->GetProjectionSurface()->F(u,v,px,py,pz);
+ patch_->GetProjectionSurface()->GetUnitNormal(u,v,nx,ny,nz);
+
+ u = patch_->UnscaleU(u);
+ v = patch_->UnscaleV(v);
+
+ d = _Interpolate(u, v).real();
+
+ x = px + d * nx;
+ y = py + d * ny;
+ z = pz + d * nz;
+ }
+ else
+ _F(u,v,x,y,z);
+}
+
+bool BlendedPatch::
+Inverse(double x,double y,double z,double &u,double &v)
+{
+ bool result;
+
+ if (_derivative) {
+ result = patch_->GetProjectionSurface()->
+ OrthoProjectionOnSurface(x,y,z,u,v);
+
+ u = patch_->UnscaleU(u);
+ v = patch_->UnscaleV(v);
+
+ double tol = 1.e-12;
+ if ((u > - tol) && (u < 1. + tol) && (v > - tol) && (v < 1. + tol))
+ result = true;
+ else
+ result = false;
+ }
+ else
+ result = _Inverse(x,y,z,u,v);
+
+ return result;
+}
+*/
+
+void BlendedPatch::
+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;
+
+ u = patch_->RescaleU(u);
+ v = patch_->RescaleV(v);
+
+ patch_->GetProjectionSurface()->F(u,v,px,py,pz);
+ patch_->GetProjectionSurface()->GetUnitNormal(u,v,nx,ny,nz);
+ patch_->GetProjectionSurface()->Dfdu(u,v,pxu,pyu,pzu);
+ patch_->GetProjectionSurface()->Dndu(u,v,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 BlendedPatch::
+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;
+
+ u = patch_->RescaleU(u);
+ v = patch_->RescaleV(v);
+
+ patch_->GetProjectionSurface()->F(u,v,px,py,pz);
+ patch_->GetProjectionSurface()->GetUnitNormal(u,v,nx,ny,nz);
+ patch_->GetProjectionSurface()->Dfdv(u,v,pxv,pyv,pzv);
+ patch_->GetProjectionSurface()->Dndv(u,v,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 BlendedPatch::
+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;
+
+ u = patch_->RescaleU(u);
+ v = patch_->RescaleV(v);
+
+ patch_->GetProjectionSurface()->F(u,v,px,py,pz);
+ patch_->GetProjectionSurface()->GetUnitNormal(u,v,nx,ny,nz);
+ patch_->GetProjectionSurface()->Dfdu(u,v,pxu,pyu,pzu);
+ patch_->GetProjectionSurface()->Dndu(u,v,nxu,nyu,nzu);
+ patch_->GetProjectionSurface()->Dfdfdudu(u,v,pxuu,pyuu,pzuu);
+ patch_->GetProjectionSurface()->Dndndudu(u,v,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 BlendedPatch::
+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;
+
+ u = patch_->RescaleU(u);
+ v = patch_->RescaleV(v);
+
+ patch_->GetProjectionSurface()->F(u,v,px,py,pz);
+ patch_->GetProjectionSurface()->GetUnitNormal(u,v,nx,ny,nz);
+ patch_->GetProjectionSurface()->Dfdv(u,v,pxv,pyv,pzv);
+ patch_->GetProjectionSurface()->Dndv(u,v,nxv,nyv,nzv);
+ patch_->GetProjectionSurface()->Dfdfdvdv(u,v,pxvv,pyvv,pzvv);
+ patch_->GetProjectionSurface()->Dndndvdv(u,v,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 BlendedPatch::
+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;
+
+ u = patch_->RescaleU(u);
+ v = patch_->RescaleV(v);
+
+ patch_->GetProjectionSurface()->F(u,v,px,py,pz);
+ patch_->GetProjectionSurface()->GetUnitNormal(u,v,nx,ny,nz);
+ patch_->GetProjectionSurface()->Dfdu(u,v,pxu,pyu,pzu);
+ patch_->GetProjectionSurface()->Dndu(u,v,nxu,nyu,nzu);
+ patch_->GetProjectionSurface()->Dfdv(u,v,pxv,pyv,pzv);
+ patch_->GetProjectionSurface()->Dndv(u,v,nxv,nyv,nzv);
+ patch_->GetProjectionSurface()->Dfdfdudv(u,v,pxuv,pyuv,pzuv);
+ patch_->GetProjectionSurface()->Dndndudv(u,v,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 BlendedPatch::GetPou(double u, double v)
+{
+ double pouU, pouV;
+
+ if (patch_->_PI->hardEdge[3])
+ pouU = OneSidedPartitionOfUnity(0.,1.,u);
+ else if (patch_->_PI->hardEdge[1])
+ pouU = 1. - OneSidedPartitionOfUnity(0.,1.,u);
+ else
+ pouU = PartitionOfUnity(u, 0., 0.3, 0.7, 1.);
+
+ if (patch_->_PI->hardEdge[0])
+ pouV = OneSidedPartitionOfUnity(0.,1.,v);
+ else if (patch_->_PI->hardEdge[2])
+ pouV = 1. - OneSidedPartitionOfUnity(0.,1.,v);
+ else
+ pouV = PartitionOfUnity(v, 0., 0.3, 0.7, 1.);
+
+ return pouU * pouV;
+}
diff --git a/contrib/FourierModel/BlendedPatch.h b/contrib/FourierModel/BlendedPatch.h
new file mode 100755
index 0000000000000000000000000000000000000000..40d0bcf6f57d50e16a59696096b52514cd90dfad
--- /dev/null
+++ b/contrib/FourierModel/BlendedPatch.h
@@ -0,0 +1,100 @@
+#ifndef _BLENDED_PATCH_H_
+#define _BLENDED_PATCH_H_
+
+#include <fftw3.h>
+#include "Message.h"
+#include "Patch.h"
+#include "BlendOperator.h"
+#include "PartitionOfUnity.h"
+
+class BlendedPatch {
+ private:
+ Patch* patch_;
+ BlendOperator* blendOp_;
+
+ int nPatches_;
+ int patchTag_;
+
+ bool IsUPeriodic, IsVPeriodic;
+
+ // bitfield telling if we also interpolate the derivative(s)
+ int _derivative;
+ // Number of Modes in reprocessed Fourier/Chebyshev series
+ int _uM, _vM;
+ // Period Information
+ double _periodU, _periodV;
+ // Limits in the series
+ int _uMLower, _uMUpper, _vMLower, _vMUpper;
+ // 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);
+ // other internal functions
+ void _ProcessSeriesCoeff();
+ // persistent fftw data (used to avoid recomputing the FFTW plans
+ // and reallocating memory every time)
+ static int _forwardSize, _backwardSize;
+ static fftw_plan _forwardPlan, _backwardPlan;
+ static fftw_complex *_forwardData, *_backwardData;
+ void _SetForwardPlan(int n);
+ void _SetBackwardPlan(int n);
+ // This routine computes the forward FFT (ignoring the last element
+ // of the data)
+ void _ForwardFft(int n, std::complex<double> *fftData);
+ // This routine computes the inverse FFT (ignoring the last element
+ // in the array), returns values in entire array (by using the
+ // periodic extension)
+ void _BackwardFft(int n, std::complex<double> *fftData);
+
+ protected:
+
+ public:
+ BlendedPatch(Patch* patch, BlendOperator* blendOp);
+ ~BlendedPatch();
+
+ double GetPou(double u, double v);
+
+ 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 Dndu
+ (double u, double v, double &x, double &y, double &z);
+
+ void Dndv
+ (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);
+
+ int GetTag
+ () { return patchTag_; }
+};
+
+#endif
diff --git a/contrib/FourierModel/CylindricalProjectionSurface.cpp b/contrib/FourierModel/CylindricalProjectionSurface.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..ba7b49d33bcd2e7af67f5e4facb6054b9eb54b78
--- /dev/null
+++ b/contrib/FourierModel/CylindricalProjectionSurface.cpp
@@ -0,0 +1,220 @@
+#include "CylindricalProjectionSurface.h"
+
+CylindricalProjectionSurface::CylindricalProjectionSurface
+(int tag) : ProjectionSurface(1.)
+{
+ SetTag(tag);
+ SetName(std::string("cylinder"));
+
+ twoPi_ = 2 * M_PI;
+
+ O_[0] = O_[1] = O_[2] = 0.;
+
+ E0_[0] = 0.; E0_[1] = 0.; E0_[2] = 1.;
+ E1_[0] = 1.; E1_[1] = 0.; E1_[2] = 0.;
+ E2_[0] = 0.; E2_[1] = 1.; E2_[2] = 0.;
+
+ scale_[0] = scale_[1] = scale_[2] = 1.;
+}
+
+CylindricalProjectionSurface::CylindricalProjectionSurface
+(int tag, double O[3], double E0[3], double E1[3], double scale[3])
+ : ProjectionSurface(1.)
+{
+ SetTag(tag);
+ SetName(std::string("cylinder"));
+
+ twoPi_ = 2 * M_PI;
+
+ O_[0] = O[0]; O_[1] = O[1]; O_[2] = O[2];
+
+ E0_[0] = E0[0]; E0_[1] = E0[1]; E0_[2] = E0[2];
+ E1_[0] = E1[0]; E1_[1] = E1[1]; E1_[2] = E1[2];
+
+ E2_[0] = E0_[1] * E1_[2] - E0_[2] * E1_[1];
+ E2_[1] = E0_[2] * E1_[0] - E0_[0] * E1_[2];
+ E2_[2] = E0_[0] * E1_[1] - E0_[1] * E1_[0];
+
+ scale_[0] = scale[0]; scale_[1] = scale[1]; scale_[2] = scale[2];
+}
+
+void CylindricalProjectionSurface::
+F(double u, double v, double &x, double &y, double &z)
+{
+ x = O_[0] + E0_[0] * scale_[0] * v;
+ y = O_[1] + E0_[1] * scale_[0] * v;
+ z = O_[2] + E0_[2] * scale_[0] * v;
+
+ x += E1_[0] * cos(twoPi_ * (u - 0.5)) +
+ E2_[0] * sin(twoPi_ * (u - 0.5));
+ y += E1_[1] * cos(twoPi_ * (u - 0.5)) +
+ E2_[1] * sin(twoPi_ * (u - 0.5));
+ z += E1_[2] * cos(twoPi_ * (u - 0.5)) +
+ E2_[2] * sin(twoPi_ * (u - 0.5));
+}
+
+bool CylindricalProjectionSurface::
+Inverse(double x, double y, double z, double &u,double &v)
+{
+ double tol =1.e-12;
+
+ double t = (x - O_[0]) * E0_[0] +
+ (y - O_[1]) * E0_[1] +
+ (z - O_[2]) * E0_[2];
+ v = t / scale_[0];
+ double n[3];
+ n[0] = x - (O_[0] + t * E0_[0]);
+ n[1] = y - (O_[1] + t * E0_[1]);
+ n[2] = z - (O_[2] + t * E0_[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]*E2_[0]+n[1]*E2_[1]+n[2]*E2_[2]),
+ (n[0]*E1_[0]+n[1]*E1_[1]+
+ n[2]*E1_[2]));
+ u /= twoPi_;
+ u += 0.5;
+
+ if ((u > - tol) && (u < 1. + tol) && (v > - tol) && (v < 1. + tol))
+ return true;
+ else
+ return false;
+}
+
+void CylindricalProjectionSurface::
+Dfdu(double u, double v, double &x, double &y, double &z)
+{
+ x = twoPi_ *
+ (- E1_[0] * sin(twoPi_ * (u - 0.5)) +
+ E2_[0] * cos(twoPi_ * (u - 0.5)));
+ y = twoPi_ *
+ (- E1_[1] * sin(twoPi_ * (u - 0.5)) +
+ E2_[1] * cos(twoPi_ * (u - 0.5)));
+ z = twoPi_ *
+ (- E1_[2] * sin(twoPi_ * (u - 0.5)) +
+ E2_[2] * cos(twoPi_ * (u - 0.5)));
+}
+
+void CylindricalProjectionSurface::
+Dfdv(double u, double v, double &x, double &y, double &z)
+{
+ x = E0_[0] * scale_[0];
+ y = E0_[1] * scale_[0];
+ z = E0_[2] * scale_[0];
+}
+
+void CylindricalProjectionSurface::
+Dfdfdudu(double u,double v, double &x, double &y, double &z)
+{
+ x = - twoPi_ * twoPi_ *
+ (E1_[0] * cos(twoPi_ * (u - 0.5)) +
+ E2_[0] * sin(twoPi_ * (u - 0.5)));
+ y = - twoPi_ * twoPi_ *
+ (E1_[1] * cos(twoPi_ * (u - 0.5)) +
+ E2_[1] * sin(twoPi_ * (u - 0.5)));
+ z = - twoPi_ * twoPi_ *
+ (E1_[2] * cos(twoPi_ * (u - 0.5)) +
+ E2_[2] * sin(twoPi_ * (u - 0.5)));
+}
+
+void CylindricalProjectionSurface::
+Dfdfdudv(double u, double v, double &x, double &y, double &z)
+{
+ x = y = z = 0.;
+}
+
+void CylindricalProjectionSurface::
+Dfdfdvdv(double u, double v, double &x, double &y, double &z)
+{
+ x = y = z = 0.;
+}
+
+void CylindricalProjectionSurface::
+Dfdfdfdududu(double u,double v,double &x,double &y,double &z)
+{
+ x = twoPi_ * twoPi_ * twoPi_ *
+ (E1_[0] * sin(twoPi_ * (u - 0.5)) -
+ E2_[0] * cos(twoPi_ * (u - 0.5)));
+ y = twoPi_ * twoPi_ * twoPi_ *
+ (E1_[1] * sin(twoPi_ * (u - 0.5)) -
+ E2_[1] * cos(twoPi_ * (u - 0.5)));
+ z = twoPi_ * twoPi_ * twoPi_ *
+ (E1_[2] * sin(twoPi_ * (u - 0.5)) +
+ E2_[2] * cos(twoPi_ * (u - 0.5)));
+}
+
+void CylindricalProjectionSurface::
+Dfdfdfdududv(double u,double v,double &x,double &y,double &z)
+{
+ x = y = z = 0.;
+}
+
+void CylindricalProjectionSurface::
+Dfdfdfdudvdv(double u,double v,double &x,double &y,double &z)
+{
+ x = y = z = 0.;
+}
+
+void CylindricalProjectionSurface::
+Dfdfdfdvdvdv(double u,double v,double &x,double &y,double &z)
+{
+ x = y = z = 0.;
+}
+
+void CylindricalProjectionSurface::
+GetNormal(double u, double v, double &x, double &y, double &z)
+{
+ x = E1_[0] * cos(twoPi_ * (u - 0.5)) +
+ E2_[0] * sin(twoPi_ * (u - 0.5));
+ y = E1_[1] * cos(twoPi_ * (u - 0.5)) +
+ E2_[1] * sin(twoPi_ * (u - 0.5));
+ z = E1_[2] * cos(twoPi_ * (u - 0.5)) +
+ E2_[2] * sin(twoPi_ * (u - 0.5));
+}
+
+void CylindricalProjectionSurface::
+Dndu(double u, double v, double &x, double &y, double &z)
+{
+ Dfdu(u,v,x,y,z);
+}
+
+void CylindricalProjectionSurface::
+Dndv(double u, double v, double &x, double &y, double &z)
+{
+ Dfdv(u,v,x,y,z);
+}
+
+void CylindricalProjectionSurface::
+Dndndudu(double u, double v, double &x, double &y, double &z)
+{
+ Dfdfdudu(u,v,x,y,z);
+}
+
+void CylindricalProjectionSurface::
+Dndndudv(double u, double v, double &x, double &y, double &z)
+{
+ Dfdfdudv(u,v,x,y,z);
+}
+
+void CylindricalProjectionSurface::
+Dndndvdv(double u, double v, double &x, double &y, double &z)
+{
+ Dfdfdvdv(u,v,x,y,z);
+}
+
+bool CylindricalProjectionSurface::
+OrthoProjectionOnSurface(double x, double y, double z, double &u,double &v)
+{
+ return Inverse(x,y,z,u,v);
+}
+
+void CylindricalProjectionSurface::
+SetParameter(int i, double x)
+{
+}
+
+double CylindricalProjectionSurface::
+GetParameter(int i)
+{
+ return 0.;
+}
diff --git a/contrib/FourierModel/CylindricalProjectionSurface.h b/contrib/FourierModel/CylindricalProjectionSurface.h
new file mode 100755
index 0000000000000000000000000000000000000000..34fc45b58197c823418389726854bf25d50a40ef
--- /dev/null
+++ b/contrib/FourierModel/CylindricalProjectionSurface.h
@@ -0,0 +1,66 @@
+#ifndef _CYLINDRICAL_PROJECTION_SURFACE_H_
+#define _CYLINDRICAL_PROJECTION_SURFACE_H_
+
+#include <cmath>
+#include "ProjectionSurface.h"
+
+class CylindricalProjectionSurface : public ProjectionSurface {
+ private:
+ double twoPi_;
+ public:
+ CylindricalProjectionSurface
+ (int tag);
+ CylindricalProjectionSurface
+ (int tag, double O[3], double E0[3], double E1[3], double scale[3]);
+
+ virtual ~CylindricalProjectionSurface
+ () {}
+
+ // Abstract methods of ProjectionSurface
+
+ 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);
+ virtual void Dfdfdfdududu
+ (double u,double v,double &x,double &y,double &z);
+ virtual void Dfdfdfdududv
+ (double u,double v,double &x,double &y,double &z);
+ virtual void Dfdfdfdudvdv
+ (double u,double v,double &x,double &y,double &z);
+ virtual void Dfdfdfdvdvdv
+ (double u,double v,double &x,double &y,double &z);
+ virtual bool OrthoProjectionOnSurface
+ (double x, double y, double z, double &u,double &v);
+ virtual void SetParameter
+ (int i, double x);
+ virtual double GetParameter
+ (int i);
+
+ // Redefinitions for CylindricalProjectionSurface
+
+ virtual void GetNormal
+ (double u, double v, double &x, double &y, double &z);
+ 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 Dndndudu
+ (double u, double v, double &x, double &y, double &z);
+ virtual void Dndndudv
+ (double u, double v, double &x, double &y, double &z);
+ virtual void Dndndvdv
+ (double u, double v, double &x, double &y, double &z);
+};
+
+#endif
diff --git a/contrib/FourierModel/Interpolator1D.cpp b/contrib/FourierModel/Interpolator1D.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..9d51584f438ab64ae0abcab27ea2442687166229
--- /dev/null
+++ b/contrib/FourierModel/Interpolator1D.cpp
@@ -0,0 +1,602 @@
+#include <vector>
+#include <math.h>
+#include <complex>
+#include "Interpolator1D.h"
+
+FftPolyInterpolator1D::FftPolyInterpolator1D(const std::vector<double> &u,
+ const std::vector< std::complex<double> > &data,
+ int refineFactor, int polyOrder, int derivative,
+ bool storeCoefs)
+ : _refineFactor(refineFactor), _polyOrder(polyOrder), _derivative(derivative),
+ _storeCoefs(storeCoefs), _uFine(0), _fineData(0), _fineDeriv(0), _fineDeriv2(0)
+{
+ _uIntervals = u.size() - 1;
+ _uFineIntervals = _refineFactor * _uIntervals;
+
+ // min/max of input mesh
+ if(u[_uIntervals] > u[0]) {
+ _uMin = u[0];
+ _uMax = u[_uIntervals];
+ }
+ else {
+ _uMin = u[_uIntervals];
+ _uMax = u[0];
+ }
+
+ // signed length of interval
+ _uLength = u[_uIntervals] - u[0];
+
+ // Sanity checks
+ if(std::abs(data[0] - data[_uIntervals]) > 1.e-10){
+ Msg::Warning("Data is not periodic: f(%g)=(%.16g,%.16g), f(%g)=(%.16g,%.16g)",
+ u[0], data[0].real(), data[0].imag(),
+ u[_uIntervals], data[_uIntervals].real(), data[_uIntervals].imag());
+ }
+ if(fabs(u[1] - u[0]) - fabs(u[2] - u[1]) > 1.e-10){
+ Msg::Fatal("Input grid is not equispaced");
+ }
+ if(_uFineIntervals < _polyOrder){
+ Msg::Fatal("Number of intervals is smaller than polynomial order (%d < %d)",
+ _uFineIntervals, _polyOrder);
+ }
+
+ // Compute the fine grid spacing (assuming the u grid is equally spaced)
+ _hUFine = (u[1] - u[0]) / _refineFactor;
+
+ // Compute the points on the fine grid
+ _uFine = new double[_uFineIntervals + 1];
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ _uFine[i] = u[0] + i * _hUFine;
+
+ // Initialize _fineData to zero
+ _fineData = new std::complex<double>[_uFineIntervals + 1];
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ _fineData[i] = 0.;
+
+ if(_storeCoefs){
+ // number of local interpolating polynomials
+ _nLocPol = 1;
+ int flag = 0;
+ while ((flag + _polyOrder) < (_uFineIntervals + 1)) {
+ _nLocPol++;
+ flag += _polyOrder;
+ }
+ if ((_uFineIntervals + 1 - flag) == 1) _nLocPol--;
+
+ // local interpolating polynomial coefficients
+ _locpol = new std::complex<double>[_nLocPol * (_polyOrder + 1)];
+ for(int i = 0; i < _nLocPol * (_polyOrder + 1); i++)
+ _locpol[i] = 0.;
+ }
+
+ std::complex<double> *forwardData = new std::complex<double>[_uIntervals+1];
+
+ if (_refineFactor == 1){ // No Fourier interpolation
+ // Copy data into _fineData
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ _fineData[i] = data[i];
+ }
+ else {
+ // Compute the forward FFT of the one-dimensional data (only uses
+ // data in i=0,...,uIntervals-1 and scales by 1/uIntervals)
+ for(int i = 0; i < _uIntervals + 1; i++)
+ forwardData[i] = data[i];
+ _ForwardFft(_uIntervals + 1, forwardData);
+
+ // Copy the Fourier coefficients into _fineData
+ int halfIndexU = _uIntervals / 2;
+ for(int i = 0; i < _uIntervals; i++){
+ int iFine = i;
+ if(i > halfIndexU)
+ iFine = _uFineIntervals + (i - _uIntervals);
+ _fineData[iFine] = forwardData[i];
+ }
+
+ // Compute inverse FFT on _fineData to interpolate the data on the
+ // fine grid
+ _BackwardFft(_uFineIntervals + 1, _fineData);
+ }
+
+ // Check if we need to interpolate the derivative(s)
+ if(_derivative){
+ if(_refineFactor == 1){
+ // No Fourier interpolation was performed, so we need to do it here...
+ for(int i = 0; i < _uIntervals + 1; i++)
+ forwardData[i] = data[i];
+ _ForwardFft(_uIntervals + 1, forwardData);
+ }
+
+ // Initialize _fineDeriv and _fineDeriv2 to zero
+ if(_derivative & 1){
+ _fineDeriv = new std::complex<double>[_uFineIntervals + 1];
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ _fineDeriv[i] = 0.;
+ if(_storeCoefs){
+ // local interpolating polynomial coefficients
+ _locpolDeriv = new std::complex<double>[_nLocPol*(_polyOrder+1)];
+ for(int i = 0; i < _nLocPol*(_polyOrder+1); i++)
+ _locpolDeriv[i] = 0.;
+ }
+ }
+ if(_derivative & 2){
+ _fineDeriv2 = new std::complex<double>[_uFineIntervals + 1];
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ _fineDeriv2[i] = 0.;
+ if(_storeCoefs){
+ // local interpolating polynomial coefficients
+ _locpolDeriv2 = new std::complex<double>[_nLocPol*(_polyOrder+1)];
+ for(int i = 0; i < _nLocPol*(_polyOrder+1); i++)
+ _locpolDeriv2[i] = 0.;
+ }
+ }
+
+ // Copy the Fourier coefficients into _fineDeriv and _fineDeriv2
+ std::complex<double> I(0., 1.);
+ int halfIndexU = _uIntervals / 2;
+ for(int i = 0; i < _uIntervals; i++){
+ int iFine = i;
+ double k = i;
+ if(i > halfIndexU){
+ k = i - _uIntervals;
+ iFine = _uFineIntervals + (i - _uIntervals);
+ }
+ // multiply by (2*pi*i*k) to get the derivative
+ if(_derivative & 1)
+ _fineDeriv[iFine] = (2. * M_PI * k * I) * forwardData[i] / _uLength;
+ // multiply by (2*pi*i*k)^2 to get the second derivative
+ if(_derivative & 2)
+ _fineDeriv2[iFine] =
+ -(4. * M_PI * M_PI * k * k) * forwardData[i] / (_uLength * _uLength);
+ }
+
+ // Perform an inverse FFT on _fineDeriv to interpolate the
+ // derivative to the fine grid
+ if(_derivative & 1)
+ _BackwardFft(_uFineIntervals + 1, _fineDeriv);
+ if(_derivative & 2)
+ _BackwardFft(_uFineIntervals + 1, _fineDeriv2);
+ }
+
+ delete [] forwardData;
+
+ if(_storeCoefs){
+ _interPolyCoeff(_locpol, _fineData, _uFineIntervals, _hUFine, _polyOrder);
+ if(_derivative & 1)
+ _interPolyCoeff(_locpolDeriv, _fineDeriv, _uFineIntervals, _hUFine,
+ _polyOrder);
+ if(_derivative & 2)
+ _interPolyCoeff(_locpolDeriv2, _fineDeriv2, _uFineIntervals, _hUFine,
+ _polyOrder);
+ }
+
+ // Initialize temp interpolation variables
+ _tmpInterp = new double[_polyOrder + 1];
+ _tmpPnts = new double[_polyOrder + 1];
+ _tmpValsReal = new double[_polyOrder + 1];
+ _tmpValsImag = new double[_polyOrder + 1];
+}
+
+FftPolyInterpolator1D::~FftPolyInterpolator1D()
+{
+ delete[] _uFine;
+ delete[] _fineData;
+ if(_storeCoefs)
+ delete[] _locpol;
+ if(_fineDeriv) {
+ delete[] _fineDeriv;
+ if(_storeCoefs)
+ delete[] _locpolDeriv;
+ }
+ if(_fineDeriv2) {
+ delete[] _fineDeriv2;
+ if(_storeCoefs)
+ delete[] _locpolDeriv2;
+ }
+
+ delete[] _tmpInterp;
+ delete[] _tmpPnts;
+ delete[] _tmpValsReal;
+ delete[] _tmpValsImag;
+}
+
+int FftPolyInterpolator1D::_forwardSize = 0;
+int FftPolyInterpolator1D::_backwardSize = 0;
+fftw_plan FftPolyInterpolator1D::_forwardPlan;
+fftw_plan FftPolyInterpolator1D::_backwardPlan;
+fftw_complex *FftPolyInterpolator1D::_forwardData = 0;
+fftw_complex *FftPolyInterpolator1D::_backwardData = 0;
+
+void FftPolyInterpolator1D::_SetForwardPlan(int n)
+{
+ if(n != _forwardSize){
+ if(_forwardSize){
+ fftw_destroy_plan(_forwardPlan);
+ fftw_free(_forwardData);
+ }
+ _forwardSize = n;
+ _forwardData = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * _forwardSize);
+ _forwardPlan = fftw_plan_dft_1d(_forwardSize, _forwardData, _forwardData,
+ FFTW_FORWARD, FFTW_ESTIMATE);
+ }
+}
+
+void FftPolyInterpolator1D::_SetBackwardPlan(int n)
+{
+ if(n != _backwardSize){
+ if(_backwardSize){
+ fftw_destroy_plan(_backwardPlan);
+ fftw_free(_backwardData);
+ }
+ _backwardSize = n;
+ _backwardData = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * _backwardSize);
+ _backwardPlan = fftw_plan_dft_1d(_backwardSize, _backwardData, _backwardData,
+ FFTW_BACKWARD, FFTW_ESTIMATE);
+ }
+}
+
+void FftPolyInterpolator1D::_ForwardFft(int n, std::complex<double> *fftData)
+{
+ // Initialize fftw plan and array (ignoring the last element of
+ // fftData, which should just be the periodic extension)
+ _SetForwardPlan(n - 1);
+ for(int i = 0; i < n - 1; i++){
+ _forwardData[i][0] = fftData[i].real();
+ _forwardData[i][1] = fftData[i].imag();
+ }
+
+ // Perform forward FFT
+ fftw_execute(_forwardPlan);
+
+ // Copy data back into fftData and scale by 1/(n - 1)
+ double s = 1. / (double)(n - 1);
+ for(int i = 0; i < n - 1; i++)
+ fftData[i] = s * std::complex<double>(_forwardData[i][0], _forwardData[i][1]);
+}
+
+void FftPolyInterpolator1D::_BackwardFft(int n, std::complex<double> *fftData)
+{
+ // Initialize fftw plan and array (ignoring last element of fftData)
+ _SetBackwardPlan(n - 1);
+ for(int i = 0; i < n - 1; i++){
+ _backwardData[i][0] = fftData[i].real();
+ _backwardData[i][1] = fftData[i].imag();
+ }
+
+ // Perform backward FFT
+ fftw_execute(_backwardPlan);
+
+ // Copy data back into fftData
+ for(int i = 0; i < n - 1; i++)
+ fftData[i] = std::complex<double>(_backwardData[i][0], _backwardData[i][1]);
+
+ // Fill in last element with copy of first element
+ fftData[n - 1] = fftData[0];
+}
+
+double FftPolyInterpolator1D::_PolyInterp(double *pnts, double *vals, double t)
+{
+ // Neville's Algorithm from Stoer and Bulirsch, Section 2.1.2
+ for(int i = 0; i <= _polyOrder; i++){
+ _tmpInterp[i] = vals[i];
+ for(int j = i-1; j >= 0; j--)
+ _tmpInterp[j] = _tmpInterp[j+1] +
+ (_tmpInterp[j+1] - _tmpInterp[j]) * (t - pnts[i]) / (pnts[i] - pnts[j]);
+ }
+ return _tmpInterp[0];
+}
+
+std::complex<double> FftPolyInterpolator1D::_Interpolate(double u, int derivative)
+{
+ if((derivative == 1 && !(_derivative & 1)) ||
+ (derivative == 2 && !(_derivative & 2)) ||
+ derivative < 0 || derivative > 2){
+ Msg::Error("Derivative data not available: check constructor call");
+ return 0.;
+ }
+
+ if(u == _uMax)
+ u = _uMin;
+
+ double epsilon = 1e-12;
+ if(u < _uMin - epsilon || u > _uMax + epsilon){
+ Msg::Error("Trying to interpolate outside interval: u=%.16g not in [%g,%g]",
+ u, _uMin, _uMax);
+ return 0.;
+ }
+
+ // Choose uIndexStart so that u is centered in the polynomial
+ // interpolation grid
+ int uIndexStart = (int)floor((u - _uFine[0]) / _hUFine) - _polyOrder / 2;
+ if(uIndexStart < 0)
+ uIndexStart = 0;
+ else if((uIndexStart + _polyOrder) > _uFineIntervals)
+ uIndexStart = _uFineIntervals - _polyOrder;
+
+ // Interpolate to find value at u
+ for(int i = 0; i <= _polyOrder; i++){
+ _tmpPnts[i] = _uFine[uIndexStart + i];
+ if(derivative == 0){
+ _tmpValsReal[i] = _fineData[uIndexStart + i].real();
+ _tmpValsImag[i] = _fineData[uIndexStart + i].imag();
+ }
+ else if(derivative == 1){
+ _tmpValsReal[i] = _fineDeriv[uIndexStart + i].real();
+ _tmpValsImag[i] = _fineDeriv[uIndexStart + i].imag();
+ }
+ else{
+ _tmpValsReal[i] = _fineDeriv2[uIndexStart + i].real();
+ _tmpValsImag[i] = _fineDeriv2[uIndexStart + i].imag();
+ }
+ }
+
+ return std::complex<double>(_PolyInterp(_tmpPnts, _tmpValsReal, u),
+ _PolyInterp(_tmpPnts, _tmpValsImag, u));
+}
+
+void polyint(double *xa, double *ya, int n, double x, double &y)
+{
+#define ind(ii) ((ii) - 1)
+ int i, m, ns = 1;
+ double dy, den, dif, dift, ho, hp, w;
+
+ double *c = new double[n];
+ double *d = new double[n];
+
+ dif = fabs(x - xa[ind(1)]);
+ for (i = 1; i <= n; i++) {
+ if ((dift = fabs(x - xa[ind(i)])) < dif) {
+ ns = i;
+ dif = dift;
+ }
+ c[ind(i)] = ya[ind(i)];
+ d[ind(i)] = ya[ind(i)];
+ }
+ y = ya[ind(ns--)];
+ for (m = 1; m < n; m++) {
+ for (i = 1; i <= n - m; i++) {
+ ho = xa[ind(i)] - x;
+ hp = xa[ind(i + m)] - x;
+ w = c[ind(i + 1)] - d[ind(i)];
+ if ((den = ho - hp) == 0) {
+ Msg::Error("Error in POLYINT");
+ return;
+ }
+ den = w / den;
+ d[ind(i)] = hp * den;
+ c[ind(i)] = ho * den;
+ }
+ if (2 * ns < (n - m))
+ dy = c[ind(ns + 1)];
+ else {
+ dy = d[ind(ns)];
+ ns--;
+ }
+ y += dy;
+ }
+ delete[] c;
+ delete[] d;
+}
+
+void polcof(double *x, double *y, int n, double *cof)
+{
+ int k;
+ double xmin;
+
+ double *xtemp = new double[n + 1];
+ double *ytemp = new double[n + 1];
+
+ for (int j = 0; j <= n; j++) {
+ xtemp[j] = x[j];
+ ytemp[j] = y[j];
+ }
+ for (int j = 0; j <= n; j++) {
+ polyint(xtemp, ytemp, n + 1 - j, 0.0, cof[j]);
+ xmin = 1.0e38;
+ k = -1;
+ for (int i = 0; i <= n - j; i++) {
+ if (fabs(xtemp[i]) < xmin) {
+ xmin = fabs(xtemp[i]);
+ k = i;
+ }
+ if (xtemp[i] != 0) ytemp[i] = (ytemp[i] - cof[j]) / xtemp[i];
+ }
+ for (int i = k + 1; i<= n - j; i++) {
+ xtemp[i - 1] = xtemp[i];
+ ytemp[i - 1] = ytemp[i];
+ }
+ }
+ delete[] xtemp;
+ delete[] ytemp;
+}
+
+void getcoeff(double *f,double *x,int n)
+{
+ double *ya = new double[n];
+ double *cof = new double[n];
+
+ for (int i = 0; i < n; i++) {
+ ya[i] = f[i];
+ }
+ polcof(x, ya, n - 1, cof);
+ for (int i = 0; i < n; i++)
+ f[i] = cof[i];
+
+ delete[] ya;
+ delete[] cof;
+}
+
+void FftPolyInterpolator1D::_interPolyCoeff(std::complex<double> *locpol,
+ std::complex<double> *u,
+ int rn, double fineh, int npoly)
+{
+ double *x = new double[npoly + 1];
+ double *freal = new double[npoly + 1];
+ double *fimag = new double[npoly + 1];
+ int rem;
+
+ int l = 0;
+ int flag = 0;
+ while ((flag + npoly) < rn) {
+ for (int i = 0; i <= npoly; i++)
+ x[i] = (flag + i) * fineh;
+ for (int i = 0; i <= npoly; i++) {
+ freal[i] = real(u[flag + i]);
+ fimag[i] = imag(u[flag + i]);
+ }
+ getcoeff(freal, x, npoly + 1);
+ getcoeff(fimag, x, npoly + 1);
+ for (int i = 0; i < npoly + 1; i++)
+ locpol[(npoly + 1) * l + i] = std::complex<double>(freal[i], fimag[i]);
+ l++;
+ flag += npoly;
+ }
+
+ rem = rn - flag;
+ if (rem > 1) {
+ for (int i = 0; i <= npoly; i++)
+ x[i] = (flag + i - (npoly + 1) + rem) * fineh;
+ for (int i=0;i<=npoly;i++) {
+ freal[i] = real(u[flag + i - (npoly + 1) + rem]);
+ fimag[i] = imag(u[flag + i - (npoly + 1) + rem]);
+ }
+ getcoeff(freal, x, npoly + 1);
+ getcoeff(fimag, x, npoly + 1);
+ for (int i = 0; i < npoly + 1; i++)
+ locpol[(npoly + 1) * l + i] = std::complex<double>(freal[i], fimag[i]);
+ }
+
+ delete[] x;
+ delete[] freal;
+ delete[] fimag;
+}
+
+std::complex<double> FftPolyInterpolator1D::_InterpolateFromCoeffs(double u,
+ int derivative)
+{
+ if((derivative == 1 && !(_derivative & 1)) ||
+ (derivative == 2 && !(_derivative & 2)) ||
+ derivative < 0 || derivative > 2){
+ Msg::Error("Derivative data not available: check constructor call");
+ return 0.;
+ }
+
+ if(u == _uMax)
+ u = _uMin;
+
+ double epsilon = 1e-12;
+ if(u < _uMin - epsilon || u > _uMax + epsilon){
+ Msg::Error("Trying to interpolate outside interval: u=%.16g not in [%g,%g]",
+ u, _uMin, _uMax);
+ return 0.;
+ }
+
+ std::complex<double> out(0.,0.);
+
+ // find the polynomial to be used for interpolating at u
+ int l = (int)floor((u - _uFine[0])/(_polyOrder * _hUFine));
+
+ if (derivative == 0) {
+ out = u * _locpol[(_polyOrder + 1) * l + _polyOrder];
+ for (int r = _polyOrder - 1; r > 0; r--)
+ out = u * (out + _locpol[(_polyOrder + 1) * l + r]);
+ out = out + _locpol[(_polyOrder + 1) * l + 0];
+ }
+ else if (derivative == 1) {
+ out = u * _locpolDeriv[(_polyOrder + 1) * l + _polyOrder];
+ for (int r = _polyOrder - 1; r > 0 ; r--)
+ out = u * (out + _locpolDeriv[(_polyOrder + 1) * l + r]);
+ out = out + _locpolDeriv[(_polyOrder + 1) * l + 0];
+ }
+ else if (derivative == 2) {
+ out = u * _locpolDeriv2[(_polyOrder + 1) * l + _polyOrder];
+ for (int r = _polyOrder - 1; r > 0; r--)
+ out = u * (out + _locpolDeriv2[(_polyOrder + 1) * l + r]);
+ out = out + _locpolDeriv2[(_polyOrder + 1) * l + 0];
+ }
+
+ return out;
+}
+
+std::complex<double> FftPolyInterpolator1D::F(double u)
+{
+ if(_storeCoefs)
+ return _InterpolateFromCoeffs(u, 0);
+ else
+ return _Interpolate(u, 0);
+}
+
+std::complex<double> FftPolyInterpolator1D::Dfdu(double u)
+{
+ if(_storeCoefs)
+ return _InterpolateFromCoeffs(u, 1);
+ else
+ return _Interpolate(u, 1);
+}
+
+std::complex<double> FftPolyInterpolator1D::Dfdfdudu(double u)
+{
+ if(_storeCoefs)
+ return _InterpolateFromCoeffs(u, 2);
+ else
+ return _Interpolate(u, 2);
+}
+
+std::complex<double> FftPolyInterpolator1D::Integrate()
+{
+ std::complex<double> value(0., 0.);
+ for(int i = 0; i < _uFineIntervals; i++){
+ value += (0.5 * _fineData[i] + 0.5 * _fineData[i + 1]) *
+ fabs(_uFine[i + 1] - _uFine[i]);
+ }
+ return value;
+}
+
+FftSplineInterpolator1D::FftSplineInterpolator1D(const std::vector<double> &u,
+ const std::vector< std::complex<double> > &data,
+ int refineFactor)
+ : FftPolyInterpolator1D(u, data, refineFactor, 3, 2)
+{
+}
+
+std::complex<double> FftSplineInterpolator1D::_SplineInterp(double *pnts,
+ std::complex<double> *vals,
+ std::complex<double> *deriv,
+ double t)
+{
+ double h = pnts[1] - pnts[0];
+ double a = (pnts[1] - t) / h;
+ double b = (t - pnts[0]) / h;
+ double a1 = a * a * a - a;
+ double b1 = b * b * b - b;
+ double h1 = (h * h) / 6.;
+ return a * vals[0] + b * vals[1] + (a1 * deriv[0] + b1 * deriv[1]) * h1;
+}
+
+std::complex<double> FftSplineInterpolator1D::F(double u)
+{
+ double pnts[2];
+ std::complex<double> vals[2], deriv[2];
+
+ double epsilon = 1e-12;
+ if(u < _uMin - epsilon || u > _uMax + epsilon){
+ Msg::Error("Trying to interpolate outside interval: u=%.16g not in [%g,%g]",
+ u, _uMin, _uMax);
+ return 0.;
+ }
+
+ // Get indices of grid points surrounding u
+ int uIndexStart = (int)floor((u - _uFine[0]) / _hUFine);
+ if(uIndexStart < 0)
+ uIndexStart = 0;
+ else if(uIndexStart > _uFineIntervals - 1)
+ uIndexStart = _uFineIntervals - 1;
+
+ pnts[0] = _uFine[uIndexStart];
+ pnts[1] = _uFine[uIndexStart + 1];
+ vals[0] = _fineData[uIndexStart];
+ vals[1] = _fineData[uIndexStart + 1];
+ deriv[0] = _fineDeriv2[uIndexStart];
+ deriv[1] = _fineDeriv2[uIndexStart + 1];
+
+ return _SplineInterp(pnts, vals, deriv, u);
+}
diff --git a/contrib/FourierModel/Interpolator1D.h b/contrib/FourierModel/Interpolator1D.h
new file mode 100755
index 0000000000000000000000000000000000000000..0db89ae1c4224aef0ae9160cdd7ba07a6c0610c9
--- /dev/null
+++ b/contrib/FourierModel/Interpolator1D.h
@@ -0,0 +1,104 @@
+#ifndef _INTERPOLATOR_1D_H_
+#define _INTERPOLATOR_1D_H_
+
+#include <vector>
+#include <complex>
+#include <fftw3.h>
+#include "Message.h"
+
+// The base class for the 1D interpolators
+class Interpolator1D {
+ public:
+ Interpolator1D(){}
+ virtual ~Interpolator1D(){}
+ // interpolates the data at u
+ virtual std::complex<double> F(double u) = 0;
+ // interpolates the first derivative of the data at u
+ virtual std::complex<double> Dfdu(double u) = 0;
+ // interpolates the second derivative of the data at u
+ virtual std::complex<double> Dfdfdudu(double u) = 0;
+};
+
+// FFT + polynomial interpolation on refined grid (assumes that the
+// input grid is equally spaced and the input data is periodic)
+class FftPolyInterpolator1D : public Interpolator1D {
+ private:
+ // refinement factor (e.g. 16; if equal to 1, we just perform
+ // standard piecewise polynomial interpolation, without any FFTs)
+ int _refineFactor;
+ // order of the interpolating polynomial
+ int _polyOrder;
+ // signed length of the interval
+ double _uLength;
+ // temporary interpolation variables
+ double *_tmpInterp, *_tmpPnts, *_tmpValsReal, *_tmpValsImag;
+ // interpolation polynomial coefficients computation routine
+ void _interPolyCoeff(std::complex<double> *locpol, std::complex<double> *u,
+ int rn,double fineh,int npoly);
+ // 1D polynomial interpolation routine
+ double _PolyInterp(double *pnts, double *vals, double t);
+ // interpolation wrappers
+ std::complex<double> _Interpolate(double u, int derivative=0);
+ std::complex<double> _InterpolateFromCoeffs(double u, int derivative=0);
+ // persistent fftw data (used to avoid recomputing the FFTW plans
+ // and reallocating memory every time)
+ static int _forwardSize, _backwardSize;
+ static fftw_plan _forwardPlan, _backwardPlan;
+ static fftw_complex *_forwardData, *_backwardData;
+ void _SetForwardPlan(int n);
+ void _SetBackwardPlan(int n);
+ protected:
+ // number of intervals in the original and in the refined mesh
+ int _uIntervals, _uFineIntervals;
+ // fine mesh spacing
+ double _hUFine;
+ // min/max of input mesh
+ double _uMin, _uMax;
+ // bitfield telling if we also interpolate the derivative(s)
+ int _derivative;
+ // flag set if we should store the polynomial coefficients
+ bool _storeCoefs;
+ // grid points of the refined mesh
+ double *_uFine;
+ // data (and its first 2 derivatives) on the refined regular grid
+ std::complex<double> *_fineData, *_fineDeriv, *_fineDeriv2;
+ // number of local interpolating polynomials
+ int _nLocPol;
+ // polynomial coefficients
+ std::complex<double> *_locpol, *_locpolDeriv, *_locpolDeriv2;
+ // This routine computes the forward FFT (ignoring the last element
+ // of the data)
+ void _ForwardFft(int n, std::complex<double> *fftData);
+ // This routine computes the inverse FFT (ignoring the last element
+ // in the array), returns values in entire array (by using the
+ // periodic extension)
+ void _BackwardFft(int n, std::complex<double> *fftData);
+ public:
+ FftPolyInterpolator1D(const std::vector<double> &u,
+ const std::vector< std::complex<double> > &data,
+ int refineFactor=16, int polyOrder=3,
+ int derivative=0, bool storeCoefs=true);
+ ~FftPolyInterpolator1D();
+ virtual std::complex<double> F(double u);
+ virtual std::complex<double> Dfdu(double u);
+ virtual std::complex<double> Dfdfdudu(double u);
+ std::complex<double> Integrate();
+};
+
+// FFT + spline interpolation on refined grid
+class FftSplineInterpolator1D : public FftPolyInterpolator1D {
+ protected:
+ // 1D spline interpolation
+ std::complex<double> _SplineInterp(double *pnts,
+ std::complex<double> *vals,
+ std::complex<double> *deriv,
+ double t);
+ public:
+ FftSplineInterpolator1D(const std::vector<double> &u,
+ const std::vector< std::complex<double> > &data,
+ int refineFactor=16);
+ ~FftSplineInterpolator1D(){}
+ virtual std::complex<double> F(double u);
+};
+
+#endif
diff --git a/contrib/FourierModel/Interpolator2D.cpp b/contrib/FourierModel/Interpolator2D.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..422708376cdea7b66341a964477a5eea8a4bc301
--- /dev/null
+++ b/contrib/FourierModel/Interpolator2D.cpp
@@ -0,0 +1,459 @@
+#include <vector>
+#include <math.h>
+#include <complex>
+#include "Interpolator2D.h"
+
+FftPolyInterpolator2D::
+FftPolyInterpolator2D(const std::vector<double> &u,
+ const std::vector<double> &v,
+ const std::vector< std::vector< std::complex<double> > > &data,
+ int refineFactor, int polyOrder, int derivative)
+ : _refineFactor(refineFactor), _polyOrder(polyOrder), _derivative(derivative),
+ _uFine(0), _vFine(0), _fineData(0), _fineDerivU(0), _fineDerivV(0),
+ _fineDerivUU(0), _fineDerivVV(0), _fineDerivUV(0)
+{
+ _uIntervals = u.size() - 1;
+ _uFineIntervals = _refineFactor * _uIntervals;
+ _vIntervals = v.size() - 1;
+ _vFineIntervals = _refineFactor * _vIntervals;
+
+ // min/max of input mesh
+ if(u[_uIntervals] > u[0]) {
+ _uMin = u[0];
+ _uMax = u[_uIntervals];
+ }
+ else {
+ _uMin = u[_uIntervals];
+ _uMax = u[0];
+ }
+ if(v[_vIntervals] > u[0]) {
+ _vMin = v[0];
+ _vMax = v[_vIntervals];
+ }
+ else {
+ _vMin = u[_vIntervals];
+ _vMax = u[0];
+ }
+
+ // signed length of interval
+ _uLength = u[_uIntervals] - u[0];
+ _vLength = v[_vIntervals] - v[0];
+
+ // Sanity checks
+ if(std::abs(data[0][0] - data[_uIntervals][0]) > 1.e-10 ||
+ std::abs(data[0][0] - data[0][_vIntervals]) > 1.e-10) {
+ Msg::Warning("Data is not periodic");
+ }
+ if(fabs(u[1] - u[0]) - fabs(u[2] - u[1]) > 1.e-10 ||
+ fabs(v[1] - v[0]) - fabs(v[2] - v[1]) > 1.e-10){
+ Msg::Fatal("Input grid is not equispaced");
+ }
+ if((_uFineIntervals < _polyOrder) || (_vFineIntervals < _polyOrder)){
+ Msg::Fatal("Number of intervals is smaller than polynomial order (%d < %d or %d < %d)",
+ _uFineIntervals, _polyOrder, _vFineIntervals, _polyOrder);
+ }
+
+ // Compute the fine grid spacing (assuming the u and v are equally spaced)
+ _hUFine = (u[1] - u[0]) / _refineFactor;
+ _hVFine = (v[1] - v[0]) / _refineFactor;
+
+ // Compute the points on the fine grid
+ _uFine = new double[_uFineIntervals + 1];
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ _uFine[i] = u[0] + i * _hUFine;
+
+ _vFine = new double[_vFineIntervals + 1];
+ for(int i = 0; i < _vFineIntervals + 1; i++)
+ _vFine[i] = v[0] + i * _hVFine;
+
+ // Initialize _fineData to zero
+ _fineData = new std::complex<double>*[_uFineIntervals + 1];
+ for(int i = 0; i < _uFineIntervals + 1; i++){
+ _fineData[i] = new std::complex<double>[_vFineIntervals + 1];
+ for(int j = 0; j < _vFineIntervals + 1; j++)
+ _fineData[i][j] = 0.;
+ }
+
+ std::complex<double> **forwardData = new std::complex<double>*[_uIntervals + 1];
+ for(int i = 0; i < _uIntervals + 1; i++)
+ forwardData[i] = new std::complex<double>[_vIntervals + 1];
+
+ if (_refineFactor == 1){ // No Fourier interpolation
+ // Copy data into _fineData
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ for(int j = 0; j < _vFineIntervals + 1; j++)
+ _fineData[i][j] = data[i][j];
+ }
+ else {
+ // Compute the forward FFT of the two dimensional data (only uses
+ // data in i=0,...,uIntervals-1 and j=0,...,vIntervals-1 and
+ // scales by 1/(uIntervals*vIntervals))
+ for(int i = 0; i < _uIntervals + 1; i++)
+ for(int j = 0; j < _vIntervals + 1; j++)
+ forwardData[i][j] = data[i][j];
+ _ForwardFft(_uIntervals + 1, _vIntervals + 1, forwardData);
+
+ // Copy the Fourier coefficients into _fineData
+ int halfIndexU = _uIntervals / 2;
+ int halfIndexV = _vIntervals / 2;
+ for(int i = 0; i < _uIntervals; i++)
+ for(int j = 0; j < _vIntervals; j++){
+ int iFine = i;
+ int jFine = j;
+ if(i > halfIndexU)
+ iFine = _uFineIntervals + (i - _uIntervals);
+ if(j > halfIndexV)
+ jFine = _vFineIntervals + (j - _vIntervals);
+ _fineData[iFine][jFine] = forwardData[i][j];
+ }
+
+ // Compute inverse FFT on _fineData to interpolate the data on the
+ // fine grid
+ _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineData);
+ }
+
+ // Check if we need to interpolate the derivative(s)
+ if(_derivative){
+ if(_refineFactor == 1){
+ // No Fourier interpolation was performed, so we need to do it here...
+ for(int i = 0; i < _uIntervals + 1; i++)
+ for(int j = 0; j < _vIntervals + 1; j++)
+ forwardData[i][j] = data[i][j];
+ _ForwardFft(_uIntervals + 1, _vIntervals + 1, forwardData);
+ }
+
+ // Initialize _fineDeriv and _fineDeriv2 to zero
+ if(_derivative & 1){
+ _fineDerivU = new std::complex<double>*[_uFineIntervals + 1];
+ _fineDerivV = new std::complex<double>*[_uFineIntervals + 1];
+ for(int i = 0; i < _uFineIntervals + 1; i++){
+ _fineDerivU[i] = new std::complex<double>[_vFineIntervals + 1];
+ _fineDerivV[i] = new std::complex<double>[_vFineIntervals + 1];
+ for(int j = 0; j < _vFineIntervals + 1; j++){
+ _fineDerivU[i][j] = 0.;
+ _fineDerivV[i][j] = 0.;
+ }
+ }
+ }
+
+ if(_derivative & 2){
+ _fineDerivUU = new std::complex<double>*[_uFineIntervals + 1];
+ _fineDerivVV = new std::complex<double>*[_uFineIntervals + 1];
+ _fineDerivUV = new std::complex<double>*[_uFineIntervals + 1];
+ for(int i = 0; i < _uFineIntervals + 1; i++){
+ _fineDerivUU[i] = new std::complex<double>[_vFineIntervals + 1];
+ _fineDerivVV[i] = new std::complex<double>[_vFineIntervals + 1];
+ _fineDerivUV[i] = new std::complex<double>[_vFineIntervals + 1];
+ for(int j = 0; j < _vFineIntervals + 1; j++){
+ _fineDerivUU[i][j] = 0.;
+ _fineDerivVV[i][j] = 0.;
+ _fineDerivUV[i][j] = 0.;
+ }
+ }
+ }
+
+ // Copy the Fourier coefficients into _fineDeriv and _fineDeriv2
+ std::complex<double> I(0., 1.);
+ int halfIndexU = _uIntervals / 2;
+ int halfIndexV = _vIntervals / 2;
+ for(int i = 0; i < _uIntervals; i++){
+ for(int j = 0; j < _vIntervals; j++){
+ int iFine = i;
+ int jFine = j;
+ double kU = i;
+ double kV = j;
+ if(i > halfIndexU){
+ kU = i - _uIntervals;
+ iFine = _uFineIntervals + (i - _uIntervals);
+ }
+ if(j > halfIndexV){
+ kV = j - _vIntervals;
+ jFine = _vFineIntervals + (j - _vIntervals);
+ }
+ if(_derivative & 1){
+ _fineDerivU[iFine][jFine] = (2. * M_PI * kU * I) * forwardData[i][j] / _uLength;
+ _fineDerivV[iFine][jFine] = (2. * M_PI * kV * I) * forwardData[i][j] / _vLength;
+ }
+ if(_derivative & 2){
+ _fineDerivUU[iFine][jFine] =
+ -(4. * M_PI * M_PI * kU * kU) * forwardData[i][j] / (_uLength * _uLength);
+ _fineDerivVV[iFine][jFine] =
+ -(4. * M_PI * M_PI * kV * kV) * forwardData[i][j] / (_vLength * _vLength);
+ _fineDerivUV[iFine][jFine] =
+ -(4. * M_PI * M_PI * kU * kV) * forwardData[i][j] / (_uLength * _vLength);
+ }
+ }
+ }
+
+ // Perform an inverse FFT on _fineDeriv to interpolate the
+ // derivative to the fine grid
+ if(_derivative & 1){
+ _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineDerivU);
+ _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineDerivV);
+ }
+ if(_derivative & 2){
+ _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineDerivUU);
+ _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineDerivVV);
+ _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineDerivUV);
+ }
+ }
+
+ for(int i = 0; i < _uIntervals + 1; i++)
+ delete [] forwardData[i];
+ delete [] forwardData;
+
+ // Initialize interpolation variables
+ _tmpSpace = new double[_polyOrder + 1];
+ _tmpValsReal = new double[_polyOrder + 1];
+ _tmpValsImag = new double[_polyOrder + 1];
+ _tmpInterpReal = new double[_polyOrder + 1];
+ _tmpInterpImag = new double[_polyOrder + 1];
+}
+
+FftPolyInterpolator2D::~FftPolyInterpolator2D()
+{
+ delete[] _uFine;
+ delete[] _vFine;
+
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ delete [] _fineData[i];
+ delete [] _fineData;
+
+ if(_fineDerivU){
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ delete [] _fineDerivU[i];
+ delete [] _fineDerivU;
+ }
+ if(_fineDerivV){
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ delete [] _fineDerivV[i];
+ delete [] _fineDerivV;
+ }
+ if(_fineDerivUU){
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ delete [] _fineDerivUU[i];
+ delete [] _fineDerivUU;
+ }
+ if(_fineDerivVV){
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ delete [] _fineDerivVV[i];
+ delete [] _fineDerivVV;
+ }
+ if(_fineDerivUV){
+ for(int i = 0; i < _uFineIntervals + 1; i++)
+ delete [] _fineDerivUV[i];
+ delete [] _fineDerivUV;
+ }
+
+ delete[] _tmpSpace;
+ delete[] _tmpValsReal;
+ delete[] _tmpValsImag;
+ delete[] _tmpInterpReal;
+ delete[] _tmpInterpImag;
+}
+
+int FftPolyInterpolator2D::_forwardSizeU = 0;
+int FftPolyInterpolator2D::_forwardSizeV = 0;
+int FftPolyInterpolator2D::_backwardSizeU = 0;
+int FftPolyInterpolator2D::_backwardSizeV = 0;
+fftw_plan FftPolyInterpolator2D::_forwardPlan;
+fftw_plan FftPolyInterpolator2D::_backwardPlan;
+fftw_complex *FftPolyInterpolator2D::_forwardData = 0;
+fftw_complex *FftPolyInterpolator2D::_backwardData = 0;
+
+void FftPolyInterpolator2D::_SetForwardPlan(int n, int m)
+{
+ if(n != _forwardSizeU || m != _forwardSizeV){
+ if(_forwardSizeU || _forwardSizeV){
+ fftw_destroy_plan(_forwardPlan);
+ fftw_free(_forwardData);
+ }
+ _forwardSizeU = n;
+ _forwardSizeV = m;
+ _forwardData = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) *
+ _forwardSizeU * _forwardSizeV);
+ _forwardPlan = fftw_plan_dft_2d(_forwardSizeU, _forwardSizeV,
+ _forwardData, _forwardData,
+ FFTW_FORWARD, FFTW_ESTIMATE);
+ }
+}
+
+void FftPolyInterpolator2D::_SetBackwardPlan(int n, int m)
+{
+ if(n != _backwardSizeU || m != _backwardSizeV){
+ if(_backwardSizeU || _backwardSizeV){
+ fftw_destroy_plan(_backwardPlan);
+ fftw_free(_backwardData);
+ }
+ _backwardSizeU = n;
+ _backwardSizeV = m;
+ _backwardData = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) *
+ _backwardSizeU * _backwardSizeV);
+ _backwardPlan = fftw_plan_dft_2d(_backwardSizeU, _backwardSizeV,
+ _backwardData, _backwardData,
+ FFTW_BACKWARD, FFTW_ESTIMATE);
+ }
+}
+
+void FftPolyInterpolator2D::_ForwardFft(int n, int m, std::complex<double> **fftData)
+{
+ // Initialize fftw plan and array (ignoring the last row and column
+ // of fftData, which should just be the periodic extension)
+ _SetForwardPlan(n - 1, m - 1);
+ int k = 0;
+ for(int i = 0; i < n - 1; i++){
+ for(int j = 0; j < m - 1; j++){
+ _forwardData[k][0] = fftData[i][j].real();
+ _forwardData[k][1] = fftData[i][j].imag();
+ k++;
+ }
+ }
+
+ // Perform forward FFT
+ fftw_execute(_forwardPlan);
+
+ // Copy data back into fftData and scale by 1/((n-1) * (m-1))
+ double s = 1. / (double)((n - 1) * (m - 1));
+ k = 0;
+ for(int i = 0; i < n - 1; i++){
+ for(int j = 0; j < m - 1; j++){
+ fftData[i][j] = s * std::complex<double>(_forwardData[k][0], _forwardData[k][1]);
+ k++;
+ }
+ }
+}
+
+void FftPolyInterpolator2D::_BackwardFft(int n, int m, std::complex<double> **fftData)
+{
+ // Initialize fftw plan and array (ignoring last row and column of fftData)
+ _SetBackwardPlan(n - 1, m - 1);
+ int k = 0;
+ for(int i = 0; i < n - 1; i++){
+ for(int j = 0; j < m - 1; j++){
+ _backwardData[k][0] = fftData[i][j].real();
+ _backwardData[k][1] = fftData[i][j].imag();
+ k++;
+ }
+ }
+
+ // Perform backward FFT
+ fftw_execute(_backwardPlan);
+
+ // Copy data back into fftData
+ k = 0;
+ for(int i = 0; i < n - 1; i++){
+ for(int j = 0; j < m - 1; j++){
+ fftData[i][j] = std::complex<double>(_backwardData[k][0], _backwardData[k][1]);
+ k++;
+ }
+ }
+
+ // Fill in last row and column with copy of first row and column
+ for(int i = 0; i <= n - 1; i++)
+ fftData[i][m - 1] = fftData[i][0];
+ for(int j = 0; j <= m - 1; j++)
+ fftData[n - 1][j] = fftData[0][j];
+}
+
+double FftPolyInterpolator2D::_PolyInterp(double start, double h, const double *vals,
+ double t, double *space)
+{
+ // Neville's Algorithm from Stoer and Bulirsch, Section 2.1.2
+ double tmp = (t - start)/h;
+ for(int i = 0; i <= _polyOrder; i++){
+ space[i] = vals[i];
+ for(int j = i - 1; j >= 0; j--)
+ space[j] = space[j + 1] + (space[j + 1] - space[j]) * ((tmp - i) / (i - j));
+ }
+ return space[0];
+}
+
+std::complex<double> FftPolyInterpolator2D::_Interpolate(double u, double v,
+ int uDer, int vDer)
+{
+ 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 < _uMin - epsilon || u > _uMax + epsilon){
+ Msg::Error("Trying to interpolate outside interval: (u,v)=(%.16g,%.16g) "
+ "not in [%g,%g]x[%g,%g]", u, v, _uMin, _uMax, _vMin, _vMax);
+ return 0.;
+ }
+
+ // Choose uIndexStart and vIndexStart so that (u,v) is centered in
+ // the polynomial interpolation grid
+ int uIndexStart = (int)floor((u - _uFine[0]) / _hUFine) - _polyOrder / 2;
+ if(uIndexStart < 0)
+ uIndexStart = 0;
+ else if((uIndexStart + _polyOrder) > _uFineIntervals)
+ uIndexStart = _uFineIntervals - _polyOrder;
+
+ int vIndexStart = (int)floor((v - _vFine[0]) / _hVFine) - _polyOrder / 2;
+ if(vIndexStart < 0)
+ vIndexStart = 0;
+ else if((vIndexStart + _polyOrder) > _vFineIntervals)
+ vIndexStart = _vFineIntervals - _polyOrder;
+
+ // Interpolate to find value at (u,v)
+ double start = _vFine[vIndexStart];
+ for(int i = 0; i <= _polyOrder; i++){
+ for(int j = 0; j <= _polyOrder; j++){
+ std::complex<double> tmp;
+ if(uDer == 0 && vDer == 0)
+ tmp = _fineData[uIndexStart + i][vIndexStart + j];
+ else if(uDer == 1 && vDer == 0)
+ tmp = _fineDerivU[uIndexStart + i][vIndexStart + j];
+ else if(uDer == 0 && vDer == 1)
+ tmp = _fineDerivV[uIndexStart + i][vIndexStart + j];
+ else if(uDer == 2 && vDer == 0)
+ tmp = _fineDerivUU[uIndexStart + i][vIndexStart + j];
+ else if(uDer == 0 && vDer == 2)
+ tmp = _fineDerivVV[uIndexStart + i][vIndexStart + j];
+ else
+ tmp = _fineDerivUV[uIndexStart + i][vIndexStart + j];
+ _tmpValsReal[j] = tmp.real();
+ _tmpValsImag[j] = tmp.imag();
+ }
+ _tmpInterpReal[i] = _PolyInterp(start, _hVFine, _tmpValsReal, v, _tmpSpace);
+ _tmpInterpImag[i] = _PolyInterp(start, _hVFine, _tmpValsImag, v, _tmpSpace);
+ }
+
+ start = _uFine[uIndexStart];
+ return std::complex<double>(_PolyInterp(start, _hUFine, _tmpInterpReal, u, _tmpSpace),
+ _PolyInterp(start, _hUFine, _tmpInterpImag, u, _tmpSpace));
+}
+
+std::complex<double> FftPolyInterpolator2D::F(double u, double v)
+{
+ return _Interpolate(u, v, 0, 0);
+}
+
+std::complex<double> FftPolyInterpolator2D::Dfdu(double u, double v)
+{
+ return _Interpolate(u, v, 1, 0);
+}
+
+std::complex<double> FftPolyInterpolator2D::Dfdv(double u, double v)
+{
+ return _Interpolate(u, v, 0, 1);
+}
+
+std::complex<double> FftPolyInterpolator2D::Dfdfdudu(double u, double v)
+{
+ return _Interpolate(u, v, 2, 0);
+}
+
+std::complex<double> FftPolyInterpolator2D::Dfdfdvdv(double u, double v)
+{
+ return _Interpolate(u, v, 0, 2);
+}
+
+std::complex<double> FftPolyInterpolator2D::Dfdfdudv(double u, double v)
+{
+ return _Interpolate(u, v, 1, 1);
+}
diff --git a/contrib/FourierModel/Interpolator2D.h b/contrib/FourierModel/Interpolator2D.h
new file mode 100755
index 0000000000000000000000000000000000000000..df44f4c7d3db85aa06230e42ccdc829d8978954d
--- /dev/null
+++ b/contrib/FourierModel/Interpolator2D.h
@@ -0,0 +1,109 @@
+#ifndef _INTERPOLATOR_2D_H_
+#define _INTERPOLATOR_2D_H_
+
+#include <vector>
+#include <complex>
+#include <fftw3.h>
+#include "Message.h"
+
+// The base class for the 2D interpolators
+class Interpolator2D {
+ public:
+ Interpolator2D(){}
+ virtual ~Interpolator2D(){}
+ // interpolates the data at u,v
+ virtual std::complex<double> F(double u, double v) = 0;
+ // interpolates the first u derivative of the data at u,v
+ virtual std::complex<double> Dfdu(double u, double v)
+ {
+ Msg::Error("First derivative not implemented for this interpolator");
+ return 0.;
+ }
+ // interpolates the first v derivative of the data at u,v
+ virtual std::complex<double> Dfdv(double u, double v)
+ {
+ Msg::Error("First derivative not implemented for this interpolator");
+ return 0.;
+ }
+ // interpolates the second u derivative of the data at u,v
+ virtual std::complex<double> Dfdfdudu(double u, double v)
+ {
+ Msg::Error("Second derivative not implemented for this interpolator");
+ return 0.;
+ }
+ // interpolates the second v derivative of the data at u,v
+ virtual std::complex<double> Dfdfdvdv(double u, double v)
+ {
+ Msg::Error("Second derivative not implemented for this interpolator");
+ return 0.;
+ }
+ // interpolates the second cross derivative of the data at u,v
+ virtual std::complex<double> Dfdfdudv(double u, double v)
+ {
+ Msg::Error("Second derivative not implemented for this interpolator");
+ return 0.;
+ }
+};
+
+// FFT + polynomial interpolation on refined grid (assumes that the
+// input grid is equally spaced and the input data is periodic)
+class FftPolyInterpolator2D : public Interpolator2D {
+ private:
+ // refinement factor (e.g. 16; if equal to 1, we just perform
+ // standard piecewise polynomial interpolation, without any FFTs)
+ int _refineFactor;
+ // order of the interpolating polynomial
+ int _polyOrder;
+ // signed length of the interval
+ double _uLength, _vLength;
+ // temporary interpolation variables
+ double *_tmpSpace, *_tmpValsReal, *_tmpValsImag, *_tmpInterpReal, *_tmpInterpImag;
+ // 1D polynomial interpolation routine
+ double _PolyInterp(double start, double h, const double *vals,
+ double t, double *space);
+ // interpolation wrapper
+ std::complex<double> _Interpolate(double u, double v, int uDer=0, int vDer=0);
+ // persistent fftw data (used to avoid recomputing the FFTW plans
+ // and reallocating memory every time)
+ static int _forwardSizeU, _forwardSizeV;
+ static int _backwardSizeU, _backwardSizeV;
+ static fftw_plan _forwardPlan, _backwardPlan;
+ static fftw_complex *_forwardData, *_backwardData;
+ void _SetForwardPlan(int n, int m);
+ void _SetBackwardPlan(int n, int m);
+ protected:
+ // number of intervals in the original and in the refined mesh
+ int _uIntervals, _vIntervals, _uFineIntervals, _vFineIntervals;
+ // fine mesh spacing
+ double _hUFine, _hVFine;
+ // min/max of input mesh
+ double _uMin, _uMax, _vMin, _vMax;
+ // bitfield telling if we also interpolate the derivative(s)
+ int _derivative;
+ // grid points of the refined mesh
+ double *_uFine, *_vFine;
+ // data (and its first 2 derivatives) on the refined regular grid
+ std::complex<double> **_fineData, **_fineDerivU, **_fineDerivV;
+ std::complex<double> **_fineDerivUU, **_fineDerivVV, **_fineDerivUV;
+ // This routine computes the forward FFT (ignoring the last row and
+ // column of the data)
+ void _ForwardFft(int n, int m, std::complex<double> **fftData);
+ // This routine computes the inverse FFT (ignoring the last row and
+ // column in the array), returns values in entire array (by using
+ // the periodic extension)
+ void _BackwardFft(int n, int m, std::complex<double> **fftData);
+ public:
+ FftPolyInterpolator2D(const std::vector<double> &u,
+ const std::vector<double> &v,
+ const std::vector< std::vector< std::complex<double> > > &data,
+ int refineFactor=16, int polyOrder=3, int derivative=0);
+ ~FftPolyInterpolator2D();
+ virtual std::complex<double> F(double u, double v);
+ virtual std::complex<double> Dfdu(double u, double v);
+ virtual std::complex<double> Dfdv(double u, double v);
+ virtual std::complex<double> Dfdfdudu(double u, double v);
+ virtual std::complex<double> Dfdfdvdv(double u, double v);
+ virtual std::complex<double> Dfdfdudv(double u, double v);
+};
+
+#endif
diff --git a/contrib/FourierModel/PartitionOfUnity.cpp b/contrib/FourierModel/PartitionOfUnity.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..13e079a5c691e3ee2215725dd3364c7fafc2e182
--- /dev/null
+++ b/contrib/FourierModel/PartitionOfUnity.cpp
@@ -0,0 +1 @@
+#include "PartitionOfUnity.h"
diff --git a/contrib/FourierModel/PartitionOfUnity.h b/contrib/FourierModel/PartitionOfUnity.h
new file mode 100755
index 0000000000000000000000000000000000000000..39f1e2e975322fa5d40b01eac95df8784c37ba37
--- /dev/null
+++ b/contrib/FourierModel/PartitionOfUnity.h
@@ -0,0 +1,136 @@
+#ifndef _PARTITION_OF_UNITY_H_
+#define _PARTITION_OF_UNITY_H_
+
+#include <cmath>
+
+inline double OneSidedPartitionOfUnity(double t0,
+ double t1,
+ double t)
+{
+ if (t <= t0) {
+ return 1.;
+ }
+ else if (t >= t1) {
+ return 0.;
+ }
+ else {
+ double x = (t - t0) / (t1 - t0);
+ return exp(2.560851702 * exp(-1. / x) / (x - 1.));
+ }
+}
+
+inline double OneSidedPartitionOfUnityDt(double t0,
+ double t1,
+ double t)
+{
+ if (t <= t0) {
+ return 0.;
+ }
+ else if (t >= t1) {
+ return 0.;
+ }
+ else {
+ // From Maple
+ return (0.2560851702e1 * pow(t - t0, -0.2e1) * (t1 - t0) *
+ exp(-0.10e1 / (t - t0) * (t1 - t0)) / ((t - t0) / (t1 - t0) - 0.10e1) -
+ 0.2560851702e1 * exp(-0.10e1 / (t - t0) * (t1 - t0)) *
+ pow((t - t0) / (t1 - t0) - 0.10e1, -0.2e1) / (t1 - t0)) *
+ exp(0.2560851702e1 * exp(-0.10e1 / (t - t0) * (t1 - t0)) /
+ ((t - t0) / (t1 - t0) - 0.10e1));
+ }
+}
+
+inline double OneSidedPartitionOfUnityDtDt(double t0,
+ double t1,
+ double t)
+{
+ if (t <= t0) {
+ return 0.;
+ }
+ else if (t >= t1) {
+ return 0.;
+ }
+ else {
+ // From Maple
+ return ((-0.5121703404E1/pow(-t0+t,3.0)*(t1-t0)*exp(-0.1E1/(-t0+t)*(t1-t0))/
+ ((-t0+t)/(t1-t0)-0.1E1)+0.2560851702E1/pow(-t0+t,4.0)*pow(t1-t0,2.0)*
+ exp(-0.1E1/(-t0+t)*(t1-t0))/((-t0+t)/(t1-t0)-0.1E1)-0.5121703404E1/pow(-t0+t,2.0)*
+ exp(-0.1E1/(-t0+t)*(t1-t0))/pow((-t0+t)/(t1-t0)-0.1E1,2.0)+0.5121703404E1*
+ exp(-0.1E1/(-t0+t)*(t1-t0))/pow((-t0+t)/(t1-t0)-0.1E1,3.0)/pow(t1-t0,2.0))*
+ exp(0.2560851702E1*exp(-0.1E1/(-t0+t)*(t1-t0))/((-t0+t)/(t1-t0)-0.1E1))+
+ pow(0.2560851702E1/pow(-t0+t,2.0)*(t1-t0)*
+ exp(-0.1E1/(-t0+t)*(t1-t0))/((-t0+t)/(t1-t0)-0.1E1)-0.2560851702E1*
+ exp(-0.1E1/(-t0+t)*(t1-t0))/pow((-t0+t)/(t1-t0)-0.1E1,2.0)/(t1-t0),2.0)*
+ exp(0.2560851702E1*exp(-0.1E1/(-t0+t)*(t1-t0))/((-t0+t)/(t1-t0)-0.1E1)));
+ }
+}
+
+// and we get the two sided POU by multiplying two One sided POUs
+inline double PartitionOfUnityInternalCall(double r,
+ double start1, double end1,
+ double end2, double start2)
+{
+ double leftPart = 1. - OneSidedPartitionOfUnity(start1, end1, r);
+ double rightPart = OneSidedPartitionOfUnity(end2, start2, r);
+ return leftPart * rightPart;
+}
+
+// and we get derivative of the the two sided POU by adding two One sided POUs
+inline double PartitionOfUnityDtInternalCall(double r,
+ double start1, double end1,
+ double end2, double start2)
+{
+ double leftPart = 1. - OneSidedPartitionOfUnity(start1, end1, r);
+ double rightPart = OneSidedPartitionOfUnity(end2, start2, r);
+ double leftPartDt = OneSidedPartitionOfUnityDt(start1, end1, r);
+ double rightPartDt = OneSidedPartitionOfUnityDt(end2, start2, r);
+ return leftPart * rightPartDt - rightPart * leftPartDt;
+}
+
+// and we get derivative of the the two sided POU by adding two One sided POUs
+inline double PartitionOfUnityDtDtInternalCall(double r,
+ double start1, double end1,
+ double end2, double start2)
+{
+ double leftPart = 1. - OneSidedPartitionOfUnity(start1, end1, r);
+ double rightPart = OneSidedPartitionOfUnity(end2, start2, r);
+ double leftPartDt = OneSidedPartitionOfUnityDt(start1, end1, r);
+ double rightPartDt = OneSidedPartitionOfUnityDt(end2, start2, r);
+ double leftPartDtDt = OneSidedPartitionOfUnityDtDt(start1, end1, r);
+ double rightPartDtDt = OneSidedPartitionOfUnityDtDt(end2, start2, r);
+
+ return leftPartDt * rightPartDt + leftPart * rightPartDtDt -
+ (rightPartDt * leftPartDt + rightPart * leftPartDtDt);
+}
+
+// Compute the value of a two sided partition of unity
+inline double PartitionOfUnity(double r,
+ double start1, double end1,
+ double end2, double start2)
+{
+ if(r < start1 || r > start2) { return 0.; }
+ if(r > end1 && r < end2) { return 1.; }
+ return PartitionOfUnityInternalCall(r, start1, end1, end2, start2);
+}
+
+// Compute the derivative of a two sided partition of unity
+inline double PartitionOfUnityDt(double r,
+ double start1, double end1,
+ double end2, double start2)
+{
+ if(r < start1 || r > start2) { return 0.; }
+ if(r > end1 && r < end2) { return 0.; }
+ return PartitionOfUnityDtInternalCall(r, start1, end1, end2, start2);
+}
+
+// Compute the second derivative of a two sided partition of unity
+inline double PartitionOfUnityDtDt(double r,
+ double start1, double end1,
+ double end2, double start2)
+{
+ if(r < start1 || r > start2) { return 0.; }
+ if(r > end1 && r < end2) { return 0.; }
+ return PartitionOfUnityDtDtInternalCall(r, start1, end1, end2, start2);
+}
+
+#endif
diff --git a/contrib/FourierModel/ProjectionSurface.cpp b/contrib/FourierModel/ProjectionSurface.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..a5db40119488cbc3fb2521532b5c725727587b2a
--- /dev/null
+++ b/contrib/FourierModel/ProjectionSurface.cpp
@@ -0,0 +1,355 @@
+#include <cmath>
+#include "ProjectionSurface.h"
+
+ProjectionSurface::ProjectionSurface
+(double uPeriod, double vPeriod)
+{
+ tag_ = -1;
+ name_ = std::string("default");
+ numParameters_ = 0;
+
+ uPeriod_ = uPeriod;
+ vPeriod_ = vPeriod;
+
+ if (uPeriod_ < 0.)
+ uPeriodic_ = false;
+ else
+ uPeriodic_ = true;
+
+ if (vPeriod_ < 0.)
+ vPeriodic_ = false;
+ else
+ vPeriodic_ = true;
+
+ for (int l=0;l<3;l++) {
+ O_[l] = E0_[l] = E1_[l] = E2_[l] = scale_[l] = 0.;
+ }
+}
+
+void ProjectionSurface::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 ProjectionSurface::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 ProjectionSurface::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 ProjectionSurface::
+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;
+}
+
+void ProjectionSurface::
+Dndndudu(double u, double v, double &x, double &y, double &z)
+{
+ double n[3],dfdu[3],dfdv[3],dfdfdudu[3],dfdfdudv[3],dfdfdvdv[3],
+ dfdfdfdududu[3],dfdfdfdududv[3],dfdfdfdudvdv[3],dfdfdfdvdvdv[3],
+ dndu[3],dndndudu[3];
+
+ GetNormal(u,v,n[0],n[1],n[2]);
+ 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]);
+ Dfdfdfdududu(u,v,dfdfdfdududu[0],dfdfdfdududu[1],dfdfdfdududu[2]);
+ Dfdfdfdududv(u,v,dfdfdfdududv[0],dfdfdfdududv[1],dfdfdfdududv[2]);
+ Dfdfdfdudvdv(u,v,dfdfdfdudvdv[0],dfdfdfdudvdv[1],dfdfdfdudvdv[2]);
+ Dfdfdfdvdvdv(u,v,dfdfdfdvdvdv[0],dfdfdfdvdvdv[1],dfdfdfdvdvdv[2]);
+
+ 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];
+
+ dndndudu[0] = dfdfdfdududu[1] * dfdv[2] + dfdfdudu[1] * dfdfdudv[2] +
+ dfdfdudu[1] * dfdfdudv[2] + dfdu[1] * dfdfdfdududv[2] -
+ dfdfdfdududu[2] * dfdv[1] - dfdfdudu[2] * dfdfdudv[1] -
+ dfdfdudu[2] * dfdfdudv[1] - dfdu[2] * dfdfdfdududv[1];
+ dndndudu[1] = dfdfdfdududu[2] * dfdv[0] + dfdfdudu[2] * dfdfdudv[0] +
+ dfdfdudu[2] * dfdfdudv[0] + dfdu[2] * dfdfdfdududv[0] -
+ dfdfdfdududu[0] * dfdv[2] - dfdfdudu[0] * dfdfdudv[2] -
+ dfdfdudu[0] * dfdfdudv[2] - dfdu[0] * dfdfdfdududv[2];
+ dndndudu[2] = dfdfdfdududu[0] * dfdv[1] + dfdfdudu[0] * dfdfdudv[1] +
+ dfdfdudu[0] * dfdfdudv[1] + dfdu[0] * dfdfdfdududv[1] -
+ dfdfdfdududu[1] * dfdv[0] - dfdfdudu[1] * dfdfdudv[0] -
+ dfdfdudu[1] * dfdfdudv[0] - dfdu[1] * dfdfdfdududv[0];
+
+ double nDotN = n[0]*n[0]+n[1]*n[1]+n[2]*n[2];
+ double nDotNu = n[0] * dndu[0] + n[1] * dndu[1] + n[2] * dndu[2];
+ double nuDotNu = dndu[0] * dndu[0] + dndu[1] * dndu[1] + dndu[2] * dndu[2];
+ double nDotNuu = n[0] * dndndudu[0] + n[1] * dndndudu[1] +
+ n[2] * dndndudu[2];
+
+ double norm = sqrt(nDotN);
+ double normCubed = nDotN * norm;
+ double normToFive = nDotN * normCubed;
+
+ x = (nDotN * dndndudu[0] - 2 * nDotNu * dndu[0] -
+ (nuDotNu + nDotNuu) * n[0]) / normCubed +
+ 3 * nDotNu * nDotNu * n[0] / normToFive;
+ y = (nDotN * dndndudu[1] - 2 * nDotNu * dndu[1] -
+ (nuDotNu + nDotNuu) * n[1]) / normCubed +
+ 3 * nDotNu * nDotNu * n[1] / normToFive;
+ z = (nDotN * dndndudu[2] - 2 * nDotNu * dndu[2] -
+ (nuDotNu + nDotNuu) * n[2]) / normCubed +
+ 3 * nDotNu * nDotNu * n[2] / normToFive;
+}
+
+void ProjectionSurface::
+Dndndudv(double u, double v, double &x, double &y, double &z)
+{
+ double n[3],dfdu[3],dfdv[3],dfdfdudu[3],dfdfdudv[3],dfdfdvdv[3],
+ dfdfdfdududu[3],dfdfdfdududv[3],dfdfdfdudvdv[3],dfdfdfdvdvdv[3],
+ dndu[3],dndv[3],dndndudv[3];
+
+ GetNormal(u,v,n[0],n[1],n[2]);
+ 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]);
+ Dfdfdfdududu(u,v,dfdfdfdududu[0],dfdfdfdududu[1],dfdfdfdududu[2]);
+ Dfdfdfdududv(u,v,dfdfdfdududv[0],dfdfdfdududv[1],dfdfdfdududv[2]);
+ Dfdfdfdudvdv(u,v,dfdfdfdudvdv[0],dfdfdfdudvdv[1],dfdfdfdudvdv[2]);
+ Dfdfdfdvdvdv(u,v,dfdfdfdvdvdv[0],dfdfdfdvdvdv[1],dfdfdfdvdvdv[2]);
+
+ 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];
+
+ 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];
+
+ dndndudv[0] = dfdfdfdududv[1] * dfdv[2] + dfdfdudu[1] * dfdfdvdv[2] +
+ dfdfdudv[1] * dfdfdudv[2] + dfdu[1] * dfdfdfdudvdv[2] -
+ dfdfdfdududv[2] * dfdv[1] - dfdfdudu[2] * dfdfdvdv[1] -
+ dfdfdudv[2] * dfdfdudv[1] - dfdu[2] * dfdfdfdudvdv[1];
+ dndndudv[1] = dfdfdfdududv[2] * dfdv[0] + dfdfdudu[2] * dfdfdvdv[0] +
+ dfdfdudv[2] * dfdfdudv[0] + dfdu[2] * dfdfdfdudvdv[0] -
+ dfdfdfdududv[0] * dfdv[2] - dfdfdudu[0] * dfdfdvdv[2] -
+ dfdfdudv[0] * dfdfdudv[2] - dfdu[0] * dfdfdfdudvdv[2];
+ dndndudv[2] = dfdfdfdududv[0] * dfdv[1] + dfdfdudu[0] * dfdfdvdv[1] +
+ dfdfdudv[0] * dfdfdudv[1] + dfdu[0] * dfdfdfdudvdv[1] -
+ dfdfdfdududv[1] * dfdv[0] - dfdfdudu[1] * dfdfdvdv[0] -
+ dfdfdudv[1] * dfdfdudv[0] - dfdu[1] * dfdfdfdudvdv[0];
+
+ double nDotN = n[0]*n[0]+n[1]*n[1]+n[2]*n[2];
+ double nDotNu = n[0] * dndu[0] + n[1] * dndu[1] + n[2] * dndu[2];
+ double nDotNv = n[0] * dndv[0] + n[1] * dndv[1] + n[2] * dndv[2];
+ double nuDotNv = dndu[0] * dndv[0] + dndu[1] * dndv[1] + dndu[2] * dndv[2];
+ double nDotNuv = n[0] * dndndudv[0] + n[1] * dndndudv[1] +
+ n[2] * dndndudv[2];
+
+ double norm = sqrt(nDotN);
+ double normCubed = nDotN * norm;
+ double normToFive = nDotN * normCubed;
+
+ x = (nDotN * dndndudv[0] - nDotNv * dndu[0] - nDotNu * dndv[0] -
+ (nuDotNv + nDotNuv) * n[0]) / normCubed +
+ 3 * nDotNu * nDotNv * n[0] / normToFive;
+ y = (nDotN * dndndudv[1] - nDotNv * dndu[1] - nDotNu * dndv[1] -
+ (nuDotNv + nDotNuv) * n[1]) / normCubed +
+ 3 * nDotNu * nDotNv * n[1] / normToFive;
+ z = (nDotN * dndndudv[2] - nDotNv * dndu[2] - nDotNu * dndv[2] -
+ (nuDotNv + nDotNuv) * n[2]) / normCubed +
+ 3 * nDotNu * nDotNv * n[2] / normToFive;
+}
+
+void ProjectionSurface::
+Dndndvdv(double u, double v, double &x, double &y, double &z)
+{
+ double n[3],dfdu[3],dfdv[3],dfdfdudu[3],dfdfdudv[3],dfdfdvdv[3],
+ dfdfdfdududu[3],dfdfdfdududv[3],dfdfdfdudvdv[3],dfdfdfdvdvdv[3],
+ dndv[3],dndndvdv[3];
+
+ GetNormal(u,v,n[0],n[1],n[2]);
+ 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]);
+ Dfdfdfdududu(u,v,dfdfdfdududu[0],dfdfdfdududu[1],dfdfdfdududu[2]);
+ Dfdfdfdududv(u,v,dfdfdfdududv[0],dfdfdfdududv[1],dfdfdfdududv[2]);
+ Dfdfdfdudvdv(u,v,dfdfdfdudvdv[0],dfdfdfdudvdv[1],dfdfdfdudvdv[2]);
+ Dfdfdfdvdvdv(u,v,dfdfdfdvdvdv[0],dfdfdfdvdvdv[1],dfdfdfdvdvdv[2]);
+
+ 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];
+
+ dndndvdv[0] = dfdfdfdudvdv[1] * dfdv[2] + dfdfdudv[1] * dfdfdvdv[2] +
+ dfdfdudv[1] * dfdfdvdv[2] + dfdu[1] * dfdfdfdvdvdv[2] -
+ dfdfdfdudvdv[2] * dfdv[1] - dfdfdudv[2] * dfdfdvdv[1] -
+ dfdfdudv[2] * dfdfdvdv[1] - dfdu[2] * dfdfdfdvdvdv[1];
+ dndndvdv[1] = dfdfdfdudvdv[2] * dfdv[0] + dfdfdudv[2] * dfdfdvdv[0] +
+ dfdfdudv[2] * dfdfdvdv[0] + dfdu[2] * dfdfdfdvdvdv[0] -
+ dfdfdfdudvdv[0] * dfdv[2] - dfdfdudv[0] * dfdfdvdv[2] -
+ dfdfdudv[0] * dfdfdvdv[2] - dfdu[0] * dfdfdfdvdvdv[2];
+ dndndvdv[2] = dfdfdfdudvdv[0] * dfdv[1] + dfdfdudv[0] * dfdfdvdv[1] +
+ dfdfdudv[0] * dfdfdvdv[1] + dfdu[0] * dfdfdfdvdvdv[1] -
+ dfdfdfdudvdv[1] * dfdv[0] - dfdfdudv[1] * dfdfdvdv[0] -
+ dfdfdudv[1] * dfdfdvdv[0] - dfdu[1] * dfdfdfdvdvdv[0];
+
+ double nDotN = n[0]*n[0]+n[1]*n[1]+n[2]*n[2];
+ double nDotNv = n[0] * dndv[0] + n[1] * dndv[1] + n[2] * dndv[2];
+ double nvDotNv = dndv[0] * dndv[0] + dndv[1] * dndv[1] + dndv[2] * dndv[2];
+ double nDotNvv = n[0] * dndndvdv[0] + n[1] * dndndvdv[1] +
+ n[2] * dndndvdv[2];
+
+ double norm = sqrt(nDotN);
+ double normCubed = nDotN * norm;
+ double normToFive = nDotN * normCubed;
+
+ x = (nDotN * dndndvdv[0] - 2 * nDotNv * dndv[0] -
+ (nvDotNv + nDotNvv) * n[0]) / normCubed +
+ 3 * nDotNv * nDotNv * n[0] / normToFive;
+ y = (nDotN * dndndvdv[1] - 2 * nDotNv * dndv[1] -
+ (nvDotNv + nDotNvv) * n[1]) / normCubed +
+ 3 * nDotNv * nDotNv * n[1] / normToFive;
+ z = (nDotN * dndndvdv[2] - 2 * nDotNv * dndv[2] -
+ (nvDotNv + nDotNvv) * n[2]) / normCubed +
+ 3 * nDotNv * nDotNv * n[2] / normToFive;
+}
+
+void ProjectionSurface::
+Rotate(double A0, double A1, double A2)
+{
+ double tol = 1.e-14;
+
+ A0 = A0 * M_PI / 180.;
+ A1 = A1 * M_PI / 180.;
+ A2 = A2 * M_PI / 180.;
+
+ double cosA0 = cos(A0);
+ double sinA0 = sin(A0);
+ double cosA1 = cos(A1);
+ double sinA1 = sin(A1);
+ double cosA2 = cos(A2);
+ double sinA2 = sin(A2);
+
+ double c[3], E0[3], E1[3], E2[3];
+
+ for (int l = 0; l < 3; l++) {
+ E0[l] = E0_[l];
+ E1[l] = E1_[l];
+ E2[l] = E2_[l];
+ }
+
+ c[0] = cosA1 * cosA2;
+ c[1] = cosA1 * sinA2;
+ c[2] = - sinA1;
+
+ for (int l = 0; l < 3; l++)
+ E0_[l] = c[0] * E0[l] + c[1] * E1[l] + c[2] * E2[l];
+
+ c[0] = sinA0 * sinA1 * cosA2 - cosA0 * sinA2;
+ c[1] = sinA0 * sinA1 * sinA2 + cosA0 * cosA2;;
+ c[2] = sinA0 * cosA1;
+
+ for (int l = 0; l < 3; l++)
+ E1_[l] = c[0] * E0[l] + c[1] * E1[l] + c[2] * E2[l];
+
+ c[0] = cosA0 * sinA1 * cosA2 + sinA0 * sinA2;
+ c[1] = cosA0 * sinA1 * sinA2 - sinA0 * cosA2;;
+ c[2] = cosA0 * cosA1;
+
+ for (int l = 0; l < 3; l++)
+ E2_[l] = c[0] * E0[l] + c[1] * E1[l] + c[2] * E2[l];
+}
+
+void ProjectionSurface::
+Translate(double D0, double D1, double D2)
+{
+ O_[0] += D0;
+ O_[1] += D1;
+ O_[2] += D2;
+}
+
+void ProjectionSurface::
+Rescale(double S0, double S1, double S2)
+{
+ scale_[0] *= S0;
+ scale_[1] *= S1;
+ scale_[2] *= S2;
+}
diff --git a/contrib/FourierModel/ProjectionSurface.h b/contrib/FourierModel/ProjectionSurface.h
new file mode 100755
index 0000000000000000000000000000000000000000..72a5573ad0976bc5bfe1a5e26d3437cd11fb5bf4
--- /dev/null
+++ b/contrib/FourierModel/ProjectionSurface.h
@@ -0,0 +1,116 @@
+#ifndef _PROJECTION_SURFACE_H_
+#define _PROJECTION_SURFACE_H_
+
+#include <string>
+
+class ProjectionSurface {
+ private:
+ int tag_;
+ std::string name_;
+ protected:
+ int numParameters_;
+ double uPeriod_, vPeriod_;
+ bool uPeriodic_, vPeriodic_;
+ double O_[3], E0_[3], E1_[3], E2_[3], scale_[3];
+ public:
+ ProjectionSurface
+ (double uPeriod = -1., double vPeriod = -1.);
+
+ virtual ~ProjectionSurface() {}
+
+ inline int GetTag() { return tag_; }
+ inline void SetTag(int tag) { tag_ = tag; }
+
+ inline std::string GetName() { return name_; }
+ inline void SetName(std::string name) { name_ = name; }
+
+ inline int GetNumParameters() { return numParameters_; }
+
+ // Public access functions
+
+ inline void GetOrigin
+ (double &x, double &y, double &z)
+ {
+ x = O_[0]; y = O_[1]; z = O_[2];
+ }
+ inline void GetE0
+ (double &x, double &y, double &z)
+ {
+ x = E0_[0]; y = E0_[1]; z = E0_[2];
+ }
+ inline void GetE1
+ (double &x, double &y, double &z)
+ {
+ x = E1_[0]; y = E1_[1]; z = E1_[2];
+ }
+ inline void GetE2
+ (double &x, double &y, double &z)
+ {
+ x = E2_[0]; y = E2_[1]; z = E2_[2];
+ }
+ inline void GetScale
+ (double &x, double &y, double &z)
+ {
+ x = scale_[0]; y = scale_[1]; z = scale_[2];
+ }
+
+ // These are the virtual functions that must be provided
+ // by all derived projection surfaces
+
+ 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;
+ virtual void Dfdfdfdududu
+ (double u,double v,double &x,double &y,double &z) = 0;
+ virtual void Dfdfdfdududv
+ (double u,double v,double &x,double &y,double &z) = 0;
+ virtual void Dfdfdfdudvdv
+ (double u,double v,double &x,double &y,double &z) = 0;
+ virtual void Dfdfdfdvdvdv
+ (double u,double v,double &x,double &y,double &z) = 0;
+ virtual bool OrthoProjectionOnSurface
+ (double x, double y, double z, double &u,double &v) = 0;
+ virtual void SetParameter
+ (int i, double x) = 0;
+ virtual double GetParameter
+ (int i) = 0;
+
+ // These functions may also be provided by the derived
+ // projection surfaces (usually for better performance),
+ // but they don't have to
+
+ virtual void GetNormal
+ (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 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 Dndndudu
+ (double u, double v, double &x, double &y, double &z);
+ virtual void Dndndudv
+ (double u, double v, double &x, double &y, double &z);
+ virtual void Dndndvdv
+ (double u, double v, double &x, double &y, double &z);
+
+ virtual void Rotate
+ (double A0, double A1, double A2);
+ virtual void Translate
+ (double D0, double D1, double D2);
+ virtual void Rescale
+ (double S0, double S1, double S2);
+};
+
+#endif
diff --git a/contrib/FourierModel/RevolvedParabolaProjectionSurface.cpp b/contrib/FourierModel/RevolvedParabolaProjectionSurface.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..1e55997506c3aafafd404d226247e4a2ad2b4dc2
--- /dev/null
+++ b/contrib/FourierModel/RevolvedParabolaProjectionSurface.cpp
@@ -0,0 +1,266 @@
+#include "RevolvedParabolaProjectionSurface.h"
+
+RevolvedParabolaProjectionSurface::RevolvedParabolaProjectionSurface
+(int tag) : ProjectionSurface(1.)
+{
+ SetTag(tag);
+ SetName(std::string("revolvedParabola"));
+
+ twoPi_ = 2 * M_PI;
+
+ R_ = 1.;
+ K_[0] = K_[1] = 0.5;
+
+ numParameters_ = 3;
+
+ O_[0] = O_[1] = O_[2] = 0.;
+
+ E0_[0] = 0.; E0_[1] = 0.; E0_[2] = 1.;
+ E1_[0] = 1.; E1_[1] = 0.; E1_[2] = 0.;
+ E2_[0] = 0.; E2_[1] = 1.; E2_[2] = 0.;
+
+ scale_[0] = scale_[1] = scale_[2] = 1.;
+}
+
+RevolvedParabolaProjectionSurface::RevolvedParabolaProjectionSurface
+(int tag, double O[3], double E0[3], double E1[3], double scale[3],
+ double R, double K[2]) : ProjectionSurface(1.)
+{
+ SetTag(tag);
+ SetName(std::string("revolvedParabola"));
+
+ twoPi_ = 2 * M_PI;
+
+ R_ = R;
+ K_[0] = K[0]; K_[1] = K[1];
+
+ O_[0] = O[0]; O_[1] = O[1]; O_[2] = O[2];
+
+ E0_[0] = E0[0]; E0_[1] = E0[1]; E0_[2] = E0[2];
+ E1_[0] = E1[0]; E1_[1] = E1[1]; E1_[2] = E1[2];
+
+ E2_[0] = E0_[1] * E1_[2] - E0_[2] * E1_[1];
+ E2_[1] = E0_[2] * E1_[0] - E0_[0] * E1_[2];
+ E2_[2] = E0_[0] * E1_[1] - E0_[1] * E1_[0];
+
+ scale_[0] = scale[0]; scale_[1] = scale[1]; scale_[2] = scale[2];
+}
+
+void RevolvedParabolaProjectionSurface::F
+(double u, double v, double &x, double &y, double &z)
+{
+ x = O_[0] + (R_ + K_[0] * (v - 0.5)) *
+ (E1_[0] * cos(twoPi_ * (u - 0.5)) + E2_[0] * sin(twoPi_ * (u - 0.5))) +
+ K_[1] * (v - 0.5) * (v - 0.5) * E0_[0];
+ y = O_[1] + (R_ + K_[0] * (v - 0.5)) *
+ (E1_[1] * cos(twoPi_ * (u - 0.5)) + E2_[1] * sin(twoPi_ * (u - 0.5))) +
+ K_[1] * (v - 0.5) * (v - 0.5) * E0_[1];
+ z = O_[2] + (R_ + K_[0] * (v - 0.5)) *
+ (E1_[2] * cos(twoPi_ * (u - 0.5)) + E2_[2] * sin(twoPi_ * (u - 0.5))) +
+ K_[1] * (v - 0.5) * (v - 0.5) * E0_[2];
+}
+
+bool RevolvedParabolaProjectionSurface::Inverse
+(double x, double y, double z, double &u,double &v)
+{
+ double R[3];
+ R[0] = x - O_[0];
+ R[1] = y - O_[1];
+ R[2] = z - O_[2];
+
+ double RdotT = 0., RdotNcT = 0.;
+ for (int i=0;i<3;i++) {
+ RdotT += R[i] * E1_[i];
+ RdotNcT += R[i] * E2_[i];
+ }
+
+ u = atan2(RdotNcT,RdotT);
+ u /= twoPi_;
+ u += 0.5;
+
+ v = sqrt(RdotT * RdotT + RdotNcT * RdotNcT);
+ v -= R_;
+ v /= K_[0];
+ v += 0.5;
+}
+
+void RevolvedParabolaProjectionSurface::Dfdu
+(double u, double v, double &x, double &y, double &z)
+{
+ x = twoPi_ * (R_ + K_[0] * (v - 0.5)) *
+ (- E1_[0] * sin(twoPi_ * (u - 0.5)) + E2_[0] * cos(twoPi_ * (u - 0.5)));
+ y = twoPi_ * (R_ + K_[0] * (v - 0.5)) *
+ (- E1_[1] * sin(twoPi_ * (u - 0.5)) + E2_[1] * cos(twoPi_ * (u - 0.5)));
+ z = twoPi_ * (R_ + K_[0] * (v - 0.5)) *
+ (- E1_[2] * sin(twoPi_ * (u - 0.5)) + E2_[2] * cos(twoPi_ * (u - 0.5)));
+}
+
+void RevolvedParabolaProjectionSurface::Dfdv
+(double u, double v, double &x, double &y, double &z)
+{
+ x = K_[0] *
+ (E1_[0] * cos(twoPi_ * (u - 0.5)) + E2_[0] * sin(twoPi_ * (u - 0.5))) +
+ 2. * K_[1] * (v - 0.5) * E0_[0];
+ y = K_[0] *
+ (E1_[1] * cos(twoPi_ * (u - 0.5)) + E2_[1] * sin(twoPi_ * (u - 0.5))) +
+ 2. * K_[1] * (v - 0.5) * E0_[1];
+ z = K_[0] *
+ (E1_[2] * cos(twoPi_ * (u - 0.5)) + E2_[2] * sin(twoPi_ * (u - 0.5))) +
+ 2. * K_[1] * (v - 0.5) * E0_[2];
+}
+
+void RevolvedParabolaProjectionSurface::Dfdfdudu
+(double u,double v, double &x, double &y, double &z)
+{
+ x = - twoPi_ * twoPi_ * (R_ + K_[0] * (v - 0.5)) *
+ (E1_[0] * cos(twoPi_ * (u - 0.5)) + E2_[0] * sin(twoPi_ * (u - 0.5)));
+ y = - twoPi_ * twoPi_ * (R_ + K_[0] * (v - 0.5)) *
+ (E1_[1] * cos(twoPi_ * (u - 0.5)) + E2_[1] * sin(twoPi_ * (u - 0.5)));
+ z = - twoPi_ * twoPi_ * (R_ + K_[0] * (v - 0.5)) *
+ (E1_[2] * cos(twoPi_ * (u - 0.5)) + E2_[2] * sin(twoPi_ * (u - 0.5)));
+}
+
+void RevolvedParabolaProjectionSurface::Dfdfdudv
+(double u, double v, double &x, double &y, double &z)
+{
+ x = K_[0] * twoPi_ *
+ (- E1_[0] * sin(twoPi_ * (u - 0.5)) + E2_[0] * cos(twoPi_ * (u - 0.5)));
+ y = K_[0] * twoPi_ *
+ (- E1_[1] * sin(twoPi_ * (u - 0.5)) + E2_[1] * cos(twoPi_ * (u - 0.5)));
+ z = K_[0] * twoPi_ *
+ (- E1_[2] * sin(twoPi_ * (u - 0.5)) + E2_[2] * cos(twoPi_ * (u - 0.5)));
+}
+
+void RevolvedParabolaProjectionSurface::Dfdfdvdv
+(double u, double v, double &x, double &y, double &z)
+{
+ x = 2. * K_[1] * E0_[0];
+ y = 2. * K_[1] * E0_[1];
+ z = 2. * K_[1] * E0_[2];
+}
+
+void RevolvedParabolaProjectionSurface::Dfdfdfdududu
+(double u,double v,double &x,double &y,double &z)
+{
+ x = twoPi_ * twoPi_ * twoPi_ * (R_ + K_[0] * (v - 0.5)) *
+ (E1_[0] * sin(twoPi_ * (u - 0.5)) - E2_[0] * cos(twoPi_ * (u - 0.5)));
+ y = twoPi_ * twoPi_ * twoPi_ * (R_ + K_[0] * (v - 0.5)) *
+ (E1_[1] * sin(twoPi_ * (u - 0.5)) - E2_[1] * cos(twoPi_ * (u - 0.5)));
+ z = twoPi_ * twoPi_ * twoPi_ * (R_ + K_[0] * (v - 0.5)) *
+ (E1_[2] * sin(twoPi_ * (u - 0.5)) + E2_[2] * cos(twoPi_ * (u - 0.5)));
+}
+
+void RevolvedParabolaProjectionSurface::Dfdfdfdududv
+(double u,double v,double &x,double &y,double &z)
+{
+ x = - twoPi_ * twoPi_ * K_[0] *
+ (E1_[0] * cos(twoPi_ * (u - 0.5)) + E2_[0] * sin(twoPi_ * (u - 0.5)));
+ y = - twoPi_ * twoPi_ * K_[0] *
+ (E1_[1] * cos(twoPi_ * (u - 0.5)) + E2_[1] * sin(twoPi_ * (u - 0.5)));
+ z = - twoPi_ * twoPi_ * K_[0] *
+ (E1_[2] * cos(twoPi_ * (u - 0.5)) + E2_[2] * sin(twoPi_ * (u - 0.5)));
+}
+
+void RevolvedParabolaProjectionSurface::Dfdfdfdudvdv
+(double u,double v,double &x,double &y,double &z)
+{
+ x = y = z = 0.;
+}
+
+void RevolvedParabolaProjectionSurface::Dfdfdfdvdvdv
+(double u,double v,double &x,double &y,double &z)
+{
+ x = y = z = 0.;
+}
+
+bool RevolvedParabolaProjectionSurface::OrthoProjectionOnSurface
+(double x, double y, double z, double& u, double& v)
+{
+ double R[3];
+ R[0] = x - O_[0];
+ R[1] = y - O_[1];
+ R[2] = z - O_[2];
+
+ double RdotT = 0., RdotNcT = 0.;
+ for (int i=0;i<3;i++) {
+ RdotT += R[i] * E1_[i];
+ RdotNcT += R[i] * E2_[i];
+ }
+
+ u = atan2(RdotNcT,RdotT);
+ u /= twoPi_;
+ u += 0.5;
+
+ double A = - R_, B = 0.;
+ for (int i=0;i<3;i++) {
+ A += R[i] * (E1_[i] * cos(twoPi_ * (u - 0.5)) +
+ E2_[i] * sin(twoPi_ * (u - 0.5)));
+ B += R[i] * E0_[i];
+ }
+
+ double a = 2 * K_[1] * K_[1];
+ double b = K_[0] * K_[0] - 2 * K_[1] * B;
+ double c = - K_[0] * A;
+
+ std::vector<double> root = SolveCubic(a,b,c);
+
+ if (root.size()) {
+ double xP,yP,zP;
+ double minDist = 1.e12;
+ for (int i=0;i<root.size();i++) {
+ F(u,root[i] + 0.5,xP,yP,zP);
+ double dist = sqrt((x-xP)*(x-xP)+(y-yP)*(y-yP)+(z-zP)*(z-zP));
+ if (dist < minDist) {
+ minDist = dist;
+ v = root[i] + 0.5;
+ }
+ }
+ double tol =1.e-4;
+ if ((u > - tol) && (u < 1. + tol) && (v > - tol) && (v < 1. + tol))
+ return true;
+ else
+ return false;
+ }
+ else
+ return false;
+}
+
+void RevolvedParabolaProjectionSurface::GetNormal
+(double u, double v, double &x, double &y, double &z)
+{
+ x = - 2. * K_[1] * (v - 0.5) *
+ (E1_[0] * cos(twoPi_ * (u - 0.5)) + E2_[0] * sin(twoPi_ * (u - 0.5))) +
+ K_[0] * E0_[0];
+ y = - 2. * K_[1] * (v - 0.5) *
+ (E1_[1] * cos(twoPi_ * (u - 0.5)) + E2_[1] * sin(twoPi_ * (u - 0.5))) +
+ K_[0] * E0_[1];
+ z = - 2. * K_[1] * (v - 0.5) *
+ (E1_[2] * cos(twoPi_ * (u - 0.5)) + E2_[2] * sin(twoPi_ * (u - 0.5))) +
+ K_[0] * E0_[2];
+}
+
+void RevolvedParabolaProjectionSurface::
+SetParameter(int i, double x)
+{
+ switch (i) {
+ case 0:
+ R_ = x;
+ case 1:
+ K_[0] = x;
+ case 2:
+ K_[1] = x;
+ }
+}
+
+double RevolvedParabolaProjectionSurface::
+GetParameter(int i)
+{
+ switch (i) {
+ case 0:
+ return R_;
+ case 1:
+ return K_[0];
+ case 2:
+ return K_[1];
+ }
+}
diff --git a/contrib/FourierModel/RevolvedParabolaProjectionSurface.h b/contrib/FourierModel/RevolvedParabolaProjectionSurface.h
new file mode 100755
index 0000000000000000000000000000000000000000..73e7a861ae02b1498ed0aa6c47e5a272d50ca733
--- /dev/null
+++ b/contrib/FourierModel/RevolvedParabolaProjectionSurface.h
@@ -0,0 +1,60 @@
+#ifndef _REVOLVED_PARABOLA_PROJECTION_SURFACE_H_
+#define _REVOLVED_PARABOLA_PROJECTION_SURFACE_H_
+
+#include <cmath>
+#include <vector>
+#include "Utils.h"
+#include "ProjectionSurface.h"
+
+class RevolvedParabolaProjectionSurface : public ProjectionSurface {
+ private:
+ double twoPi_;
+ double R_, K_[2];
+ public:
+ RevolvedParabolaProjectionSurface
+ (int tag);
+ RevolvedParabolaProjectionSurface
+ (int tag, double O[3], double E0[3], double E1[3], double scale[3],
+ double R, double K[2]);
+
+ virtual ~RevolvedParabolaProjectionSurface
+ () {}
+
+ // Abstract methods of ProjectionSurface
+
+ 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);
+ virtual void Dfdfdfdududu
+ (double u,double v,double &x,double &y,double &z);
+ virtual void Dfdfdfdududv
+ (double u,double v,double &x,double &y,double &z);
+ virtual void Dfdfdfdudvdv
+ (double u,double v,double &x,double &y,double &z);
+ virtual void Dfdfdfdvdvdv
+ (double u,double v,double &x,double &y,double &z);
+ virtual bool OrthoProjectionOnSurface
+ (double x, double y, double z, double &u,double &v);
+ virtual void SetParameter
+ (int i, double x);
+ virtual double GetParameter
+ (int i);
+
+ // Redefinitions for RevolvedParabolaProjectionSurface
+
+ virtual void GetNormal
+ (double u, double v, double &x, double &y, double &z);
+};
+
+#endif
diff --git a/contrib/FourierModel/Utils.cpp b/contrib/FourierModel/Utils.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..956c8f8d9e1927d829b1056a5e755da0856ecf57
--- /dev/null
+++ b/contrib/FourierModel/Utils.cpp
@@ -0,0 +1,87 @@
+#include "Utils.h"
+
+std::vector<double> SolveCubic(double a, double b, double c)
+{
+ std::vector<double> root;
+ // find real roots of polynomial a*x^3+ b*x +c=0
+
+ double tol=pow(10.0, -12.);
+
+ double error=1.;
+ int maxiter=100; int iter=0;
+
+ std::vector<double> initial_guess(5);
+ initial_guess[0]=0.;
+ initial_guess[1]=1.;
+ initial_guess[2]=-1.;
+ initial_guess[3]=0.5;
+ initial_guess[4]=-0.5;
+
+ bool find=0; int init=0;
+ double x0=0., x1, f, df;
+
+ while (find==0 && init< 5){
+
+ init=init+1;
+ x0=initial_guess[init-1];
+ error=1.; iter=0;
+
+ while (error>tol && iter<maxiter){
+
+ f=a*x0*x0*x0+b*x0+c;
+ df=a*3.*x0*x0+b;
+
+ x1=x0-f/df;
+
+ error=std::abs(x1-x0);
+ x0=x1;
+ iter=iter+1;
+ }
+
+ if (iter<maxiter){
+ find=1;
+ }
+ }
+
+ if (find) {
+ root.push_back(x0);
+
+ //%polynomial a*x^3+ b*x +c=(x-x0)*(a*x^2+bb*x+cc)
+ //%where bb=a*x0, cc=-c/x0 (=b-a*x^2);
+
+ double aa=a, bb=a*x0, cc=b-a*x0*x0;
+ double D=bb*bb-4.*aa*cc;
+
+ if (D>=0.){
+ root.push_back((-bb+sqrt(D))/2./aa);
+ root.push_back((-bb-sqrt(D))/2./aa);
+ }
+
+ // it seems that the solution obtained by quadratic formula
+ //is not so accurate --> used as initial guess.
+
+ if (root.size() == 3){
+ for (int i=0 ; i < 2 ; i++){
+ x0=root[i+1];
+ error=1.; iter=0;
+ while (error>tol && iter<maxiter){
+
+ f=a*x0*x0*x0+b*x0+c;
+ df=a*3.*x0*x0+b;
+
+ x1=x0-f/df;
+
+ error=std::abs(x1-x0);
+ x0=x1;
+ iter=iter+1;
+ }
+ root[i+1]=x0;
+ }
+ }
+ }
+
+ // cout << a << ' ' << b << ' ' << c << ' ' <<
+ // num << ' ' << root[0] << ' ' <<root[1] << ' ' << root[2] << endl;
+
+ return root;
+}
diff --git a/contrib/FourierModel/Utils.h b/contrib/FourierModel/Utils.h
new file mode 100755
index 0000000000000000000000000000000000000000..9a51797d7e5f154d8dd84d19dc180235ac97974a
--- /dev/null
+++ b/contrib/FourierModel/Utils.h
@@ -0,0 +1,9 @@
+#ifndef _UTILS_H_
+#define _UTILS_H_
+
+#include <cmath>
+#include <vector>
+
+std::vector<double> SolveCubic(double a, double b, double c);
+
+#endif