more changes

Fltk/GUI_Projection.cpp

@@ -854,6 +854,28 @@ void compute_cb(Fl_Widget *w, void *data)
     }
   }
 
+  // IO Test Code
+
+  char *filename = "patches.fm";
+
+  FILE *fp = fopen(filename, "w");
+  if(!fp){
+    printf("Unable to open file '%s'\n", filename);
+    return;
+  }
+
+  std::set<GFace*, GEntityLessThan>::iterator fiter;
+  for (fiter = GMODEL->firstFace(); fiter != GMODEL->lastFace(); fiter++) {
+    if ((*fiter)->getNativeType() == GEntity::FourierModel) {
+      FFace* ff = (FFace*) (*fiter);
+      ff->GetFMFace()->GetPatch()->Export(fp);
+    }
+  }
+
+  FM_Reader* reader = new FM_Reader(filename);
+
+  // End Test
+
   Draw();
 }
 
@@ -911,10 +933,15 @@ void mesh_parameterize_cb(Fl_Widget* w, void* data)
   if(!editor){
     std::vector<FProjectionFace*> faces;
     int tag = GMODEL->numFace();
+    faces.push_back(new FProjectionFace(GMODEL, ++tag,
+					new PlaneProjectionSurface(tag)));
+    faces.push_back(new FProjectionFace(GMODEL, ++tag,
+					new ParaboloidProjectionSurface(tag)));
     faces.push_back(new FProjectionFace(GMODEL, ++tag, 
 					new CylindricalProjectionSurface(tag)));
     faces.push_back(new FProjectionFace(GMODEL, ++tag,
 					new RevolvedParabolaProjectionSurface(tag)));
+
     editor = new projectionEditor(faces);
 
     for(unsigned int i = 0; i < faces.size(); i++){

Fltk/GUI_Projection.h

@@ -24,6 +24,7 @@
 #include "FM_Vertex.h"
 #include "FM_Edge.h"
 #include "FM_Face.h"
+#include "FM_Reader.h"
 
 void select_cb(Fl_Widget *w, void *data);
 void browse_cb(Fl_Widget *w, void *data);

Geo/FFace.h

@@ -42,6 +42,7 @@ class FFace : public GFace {
   ModelType getNativeType() const { return FourierModel; }
   void * getNativePtr() const {throw;} 
   virtual bool surfPeriodic(int dim) const;
+  FM_Face* GetFMFace() { return face; }
 };
 
 #endif

Geo/FProjectionFace.h

@@ -7,6 +7,8 @@
 #if defined(HAVE_FOURIER_MODEL)
 
 #include "ProjectionSurface.h"
+#include "PlaneProjectionSurface.h"
+#include "ParaboloidProjectionSurface.h"
 #include "CylindricalProjectionSurface.h"
 #include "RevolvedParabolaProjectionSurface.h"
 

Geo/GModelIO_F.cpp

@@ -19,6 +19,7 @@ void F_Internals::loadF(const char *fn)
 
 void F_Internals::buildGModel(GModel *model)
 {
+/*
   // building geom vertices
   int nVertices = _reader->GetNumVertices();
   for (int i=0;i<nVertices;i++) {
@@ -48,6 +49,7 @@ void F_Internals::buildGModel(GModel *model)
     FFace *f = new FFace(model, face, i, l_edges);
     model->add(f);
   }
+*/
 }
 
 int GModel::readF(const std::string &fn)

contrib/FourierModel/ContinuationPatch.cpp

@@ -3,13 +3,14 @@
 #include "ContinuationPatch.h"
 
 ContinuationPatch::ContinuationPatch
-(PatchInfo* PI, ProjectionSurface* ps, int derivative) 
+(PatchInfo* PI, ProjectionSurface* ps)
   : Patch(),_coeffOriginalData(0),_coeffData(0),_coeffDerivU(0),
     _coeffDerivV(0),_coeffDerivUU(0),_coeffDerivVV(0),_coeffDerivUV(0)
 {
   _PI = PI;
   _tag = _PI->tag;
-  _derivative = derivative;
+  _derivative = _PI->derivative;
+  _recompute = _PI->recompute;
 
   SetProjectionSurface(ps);
 
@@ -18,14 +19,23 @@ ContinuationPatch::ContinuationPatch
   _vMin = _PI->vMin;
   _vMax = _PI->vMax;
 
-  _periodicityU = _PI->periodic[0];
-  _periodicityV = _PI->periodic[1];
+  if (_ps->IsUPeriodic())
+    _periodicityU = 1;
+  if (_ps->IsVPeriodic())
+    _periodicityV = 1;
 
   for (int i = 0; i < 4; i++)
     _hardEdge[i] = _PI->hardEdge[i];
 
-  _periodU = (1-_PI->periodic[0]) + 1;
-  _periodV = (1-_PI->periodic[1]) + 1;
+  if (_ps->IsUPeriodic())
+    _periodU = 1.;
+  else
+    _periodU = 2.;
+
+  if (_ps->IsVPeriodic())
+    _periodV = 1.;
+  else
+    _periodV = 2.;
 
   _uModes = _PI->nModes[0];
   _vModes = _PI->nModes[1];
@@ -79,46 +89,43 @@ ContinuationPatch::ContinuationPatch
     _vMUpper = _vM;    
   }
 
+  // Initialize interpolation variables
+  _tmpCoeff = std::vector< std::complex<double> >(_vM);
+  _tmpInterp = std::vector< std::complex<double> >(_uM);
+
   // Initialize Data
   _coeffOriginalData = new std::complex<double>*[_uModes];
   for(int j = 0; j < _uModes; j++){
     _coeffOriginalData[j] = new std::complex<double>[_vModes];
     for(int k = 0; k < _vModes; k++)
-      _coeffOriginalData[j][k] = _PI->coeff[j][k];
+      _coeffOriginalData[j][k] = _PI->coeffFourier[j][k];
   }
 
   // Initialize interpolation variables
   _tmpOrigCoeff = std::vector< std::complex<double> >(_vModes);
   _tmpOrigInterp = std::vector< std::complex<double> >(_uModes);
 
-  if (_derivative)
-    _ReprocessSeriesCoeff();
-
-  /*
-  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) {
+    // Initialize Chebyshev coefficients
+    _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.;
+    }
+    // Initialize Deriv and Deriv2 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.;
-        }
+	_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];
@@ -134,214 +141,140 @@ ContinuationPatch::ContinuationPatch
         }
       }
     }
-
-    // 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.;
+    if (_recompute) {
+      _ReprocessSeriesCoeff();
+      // 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 (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];
+	    }
 	  }
-        }
+	}
       }
-    }
-    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.;
+      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.;
+	    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.;
+      
+      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];
 	    }
-	    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];
+	    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 {
-	      _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 == _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];
+	      }
 	    }
 	  }
 	}
       }
-    }
-    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.;
+      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.;
+	    if (IsUPeriodic() && !IsVPeriodic()) {
+	      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++)
-      printf("%g ",_coeffDerivVV[j][k].real());
-    printf("\n\n");
-  }
-  printf("\n");
-  */
-
-  // Initialize interpolation variables
-  _tmpCoeff = std::vector< std::complex<double> >(_vM);
-  _tmpInterp = std::vector< std::complex<double> >(_uM);
-
-  /*
-  int nU = 64;
-  int nV = 64;
-
-  double hU = 1. / (double)(nU-1);
-  double hV = 1. / (double)(nV-1);
-
-  double f_error = 0.;
-  double dfdu_error = 0.;
-  double dfdv_error = 0.;
-  double dfdfdudu_error = 0.;
-  double dfdfdudv_error = 0.;
-  double dfdfdvdv_error = 0.;
-
-  for (int j = 0; j < nU; j++)
-    for (int k = 0; k < nV; k++) {
-      double u = j * hU;
-      double v = k * hV;
-      //printf("%d %d\n",j,k);
-      std::complex<double> f = _Interpolate(u,v,0,0);
-      //printf("%d %d\n",j,k);
-      //std::complex<double> ef = 2. * v * v - 1.;
-      //std::complex<double> ef = 2. * v * v * cos(2 * M_PI * u);
-      std::complex<double> ef = cos(v) * cos(2 * M_PI * u);
-      f_error = std::max(f_error,std::abs(f - ef));
-
-      std::complex<double> dfdu = _Interpolate(u,v,1,0);
-      //std::complex<double> edfdu = 0.;
-      //std::complex<double> edfdu = - 4. * M_PI * v * v * sin(2 * M_PI * u);
-      std::complex<double> edfdu = - 2. * M_PI * cos(v) * sin(2 * M_PI * u);
-      dfdu_error = std::max(dfdu_error,std::abs(dfdu - edfdu));
-
-      std::complex<double> dfdv = _Interpolate(u,v,0,1);
-      //std::complex<double> edfdv = 4. * v;
-      //std::complex<double> edfdv = 4. * v * cos(2 * M_PI * u);
-      std::complex<double> edfdv = - sin(v) * cos(2 * M_PI * u);
-      dfdv_error = std::max(dfdv_error,std::abs(dfdv - edfdv));
-      //printf("%d %d : %g %g :: %g %g\n",j,k,dfdv.real(),edfdv.real(),
-      //   dfdv.imag(),edfdv.imag());
-
-      std::complex<double> dfdfdudu = _Interpolate(u,v,2,0);
-      //std::complex<double> edfdfdudu = 0.;
-      //std::complex<double> edfdfdudu = 
-      //- 8. * M_PI * M_PI * v * v * cos(2 * M_PI * u);
-      std::complex<double> edfdfdudu = 
-	- 4. * M_PI * M_PI * cos(v) * cos(2 * M_PI * u);
-      dfdfdudu_error = std::max(dfdfdudu_error,std::abs(dfdfdudu - edfdfdudu));
-
-      std::complex<double> dfdfdvdv = _Interpolate(u,v,0,2);
-      //std::complex<double> edfdfdvdv = 4.;
-      //std::complex<double> edfdfdvdv = 4. * cos(2 * M_PI * u);
-      std::complex<double> edfdfdvdv = - cos(v) * cos(2 * M_PI * u);
-      dfdfdvdv_error = std::max(dfdfdvdv_error,std::abs(dfdfdvdv - edfdfdvdv));
-      //printf("%d %d : %g %g :: %g %g\n",j,k,dfdfdvdv.real(),edfdfdvdv.real(),
-      //     dfdfdvdv.imag(),edfdfdvdv.imag());
-
-      std::complex<double> dfdfdudv = _Interpolate(u,v,1,1);
-      //std::complex<double> edfdfdudv = 0.;
-      //std::complex<double> edfdfdudv = - 8. * M_PI * v * sin(2 * M_PI * u);
-      std::complex<double> edfdfdudv = 2. * M_PI * sin(v) * sin(2 * M_PI * u);
-      dfdfdudv_error = std::max(dfdfdudv_error,std::abs(dfdfdudv - edfdfdudv));
+    else {
+      for (int j = 0; j < _uM; j++) {
+	for (int k = 0; k < _vModes; k++) {
+	  _coeffData[j][k] = _PI->coeffCheby[j][k];
+	  _coeffDerivU[j][k] = _PI->coeffDerivU[j][k];
+	  _coeffDerivV[j][k] = _PI->coeffDerivV[j][k];
+	  _coeffDerivUU[j][k] = _PI->coeffDerivUU[j][k];
+	  _coeffDerivUV[j][k] = _PI->coeffDerivUV[j][k];
+	  _coeffDerivVV[j][k] = _PI->coeffDerivVV[j][k];
+	}
+      }
     }
-
-  printf("F_Error = %g\n",f_error);
-  printf("Dfdu_Error = %g\n",dfdu_error);
-  printf("Dfdv_Error = %g\n",dfdv_error);
-  printf("Dfdfdudu_Error = %g\n",dfdfdudu_error);
-  printf("Dfdfdudv_Error = %g\n",dfdfdudv_error);
-  printf("Dfdfdvdv_Error = %g\n\n",dfdfdvdv_error);
-  */
+  }
 }
 
 ContinuationPatch::~ContinuationPatch()
@@ -463,12 +396,6 @@ void ContinuationPatch::_BackwardFft(int n, std::complex<double> *fftData)
 
 void ContinuationPatch::_ReprocessSeriesCoeff()
 {
-  _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() && IsVPeriodic()) {
     int uShift = (_uM-_uModes)%2 == 0 ? (_uM-_uModes)/2 : (_uM-_uModes)/2 + 1;
     int vShift = (_vM-_vModes)%2 == 0 ? (_vM-_vModes)/2 : (_vM-_vModes)/2 + 1;
@@ -713,10 +640,7 @@ F(double u, double v, double &x, double &y, double &z)
   _ps->F(u,v,px,py,pz);
   _ps->GetUnitNormal(u,v,nx,ny,nz);
 
-  if (_derivative)
-    d = _Interpolate(u, v, 0, 0).real();
-  else
-    d = _Interpolate(u, v).real();
+  d = _Interpolate(u, v).real();
 
   x = px + d * nx;
   y = py + d * ny;
@@ -886,3 +810,76 @@ double  ContinuationPatch::GetPou(double u, double v)
 
   return pouU * pouV;
 }
+
+void ContinuationPatch::Export(FILE *fp)
+{
+  double x,y,z;
+
+  fprintf(fp, "%s\n", _ps->GetName().c_str());
+  fprintf(fp, "%d\n", _ps->GetTag());
+  _ps->GetOrigin(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  _ps->GetE0(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  _ps->GetE1(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  _ps->GetScale(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  fprintf(fp, "%d\n", _ps->GetNumParameters());
+  for (int i = 0; i < _ps->GetNumParameters(); i++)
+    fprintf(fp, "%g\n",_ps->GetParameter(i));
+
+  fprintf(fp, "%s\n", "ContinuationPatch");
+  fprintf(fp, "%d\n", _tag);
+  fprintf(fp, "%d\n", _derivative);
+  fprintf(fp, "%g %g\n", _uMin, _uMax);
+  fprintf(fp, "%g %g\n", _vMin, _vMax);
+  fprintf(fp, "%d %d %d %d\n", _hardEdge[0], _hardEdge[1], _hardEdge[2],
+	  _hardEdge[3]);
+  fprintf(fp, "%d %d\n", _uModes, _vModes);
+  for (int j = 0; j < _uModes; j++) {
+    for (int k = 0; k < _vModes; k++) {
+      fprintf(fp, "%g %g\n", _coeffOriginalData[j][k].real(),
+	      _coeffOriginalData[j][k].imag());
+    }
+  }
+  fprintf(fp, "%d\n", 0);
+  fprintf(fp, "%d %d\n", _uM, _vM);
+  for (int j = 0; j < _uM; j++) {
+    for (int k = 0; k < _vM; k++) {
+      fprintf(fp, "%g %g\n", _coeffData[j][k].real(), _coeffData[j][k].imag());
+    }
+  }
+  if (_derivative) {
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffDerivU[j][k].real(), 
+		_coeffDerivU[j][k].imag());
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffDerivV[j][k].real(), 
+		_coeffDerivV[j][k].imag());
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffDerivUU[j][k].real(), 
+		_coeffDerivUU[j][k].imag());
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffDerivUV[j][k].real(), 
+		_coeffDerivUV[j][k].imag());
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffDerivVV[j][k].real(), 
+		_coeffDerivVV[j][k].imag());
+      }
+    }
+  }
+}

contrib/FourierModel/ContinuationPatch.h

@@ -9,6 +9,8 @@
 // The base class for the patches
 class ContinuationPatch : public Patch {
  protected:
+  // Do we recompute
+  int _recompute;
   // Number of Fourier Modes
   int _uModes, _vModes;
   // Number of Modes in reprocessed Fourier/Chebyshev series
@@ -50,13 +52,14 @@ class ContinuationPatch : public Patch {
   void _BackwardFft(int n, std::complex<double> *fftData);
  public:
   ContinuationPatch
-    (PatchInfo* PI, ProjectionSurface* ps, int derivative = 0);
+    (PatchInfo* PI, ProjectionSurface* ps);
   virtual ~ContinuationPatch();
 
   PatchInfo* _PI;
 
   // Abstract functions of Patch
 
+  virtual void Export(FILE *fp);
   virtual double GetPou(double u, double v);
   virtual void F(double u, double v, double &x, double &y, double &z);
   virtual bool Inverse(double x,double y,double z,double &u,double &v);

contrib/FourierModel/CylindricalProjectionSurface.cpp

@@ -43,6 +43,31 @@ CylindricalProjectionSurface::CylindricalProjectionSurface
   scale_[0] = scale[0]; scale_[1] = scale[1]; scale_[2] = scale[2];
 }
 
+CylindricalProjectionSurface::CylindricalProjectionSurface
+(int tag, double O[3], double E0[3], double E1[3], double scale[3],
+ double R, double Z)
+  : 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];
+
+  R_ = R;
+  Z_ = Z;
+}
+
 void CylindricalProjectionSurface::
 F(double u, double v, double &x, double &y, double &z)
 {

contrib/FourierModel/CylindricalProjectionSurface.h

@@ -13,6 +13,9 @@ class CylindricalProjectionSurface : public ProjectionSurface {
     (int tag);
   CylindricalProjectionSurface
     (int tag, double O[3], double E0[3], double E1[3], double scale[3]);
+  CylindricalProjectionSurface
+    (int tag, double O[3], double E0[3], double E1[3], double scale[3],
+     double R, double Z);
   CylindricalProjectionSurface(CylindricalProjectionSurface *ps) 
     : ProjectionSurface(ps) 
   {

contrib/FourierModel/ExactPatch.cpp

 #include "Message.h"
 #include "ExactPatch.h"
 
-ExactPatch::ExactPatch(PatchInfo *PI, ProjectionSurface* ps) : Patch()
+ExactPatch::ExactPatch(ProjectionSurface* ps) : Patch()
 {
-  _PI = PI;
   SetProjectionSurface(ps);
-
-  if (!strcmp(_PI->type,"plane")) {
-    _PI->periodic[0] = 0;
-    _PI->periodic[1] = 0;
-  }
-  else if (!strcmp(_PI->type,"disc")) {
-    _PI->periodic[0] = 0;
-    _PI->periodic[1] = 1;
-  }
-  else {
-    Msg::Info("Unknown exact patch type");
-    _PI->periodic[0] = 0;
-    _PI->periodic[1] = 0;
-  }
 }
 
 void ExactPatch::F(double u, double v, double &x, double &y, double &z)
@@ -60,3 +45,30 @@ double ExactPatch::GetPou(double u, double v)
 {
   return 1.;
 }
+
+void ExactPatch::Export(FILE *fp)
+{
+  double x,y,z;
+
+  fprintf(fp, "%s\n", _ps->GetName().c_str());
+  fprintf(fp, "%d\n", _ps->GetTag());
+  _ps->GetOrigin(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  _ps->GetE0(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  _ps->GetE1(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  _ps->GetScale(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  fprintf(fp, "%d\n", _ps->GetNumParameters());
+  for (int i = 0; i < _ps->GetNumParameters(); i++)
+    fprintf(fp, "%g\n",_ps->GetParameter(i));
+
+  fprintf(fp, "%s\n", "ExactPatch");
+  fprintf(fp, "%d\n", _tag);
+  fprintf(fp, "%d\n", _derivative);
+  fprintf(fp, "%g %g\n", _uMin, _uMax);
+  fprintf(fp, "%g %g\n", _vMin, _vMax);
+  fprintf(fp, "%d %d %d %d\n", _hardEdge[0], _hardEdge[1], _hardEdge[2],
+	  _hardEdge[3]);
+}

contrib/FourierModel/ExactPatch.h

@@ -8,17 +8,17 @@
 // The base class for the patches
 class ExactPatch : public Patch {
  public:
-  ExactPatch(PatchInfo* PI, ProjectionSurface* ps);
+  ExactPatch(ProjectionSurface* ps);
   virtual ~ExactPatch() {}
 
-  PatchInfo* _PI;
-
   // These are the virtual functions that must be provided by all
   // derived patches: GetPou() returns the original smooth
   // (non-normalized) cutoff on the patch; F() and Inverse() implement
   // the mapping f: (u,v)->(x,y,z) and its inverse; and the Df*() and Dn*()
   // functions return the derivatives of the mapping f and unit normal n 
   // with respect to u and v
+
+  virtual void Export(FILE *fp);
   virtual double GetPou(double u, double v);
   virtual void F(double u, double v, double &x, double &y, double &z);
   virtual bool Inverse(double x,double y,double z,double &u,double &v);

contrib/FourierModel/FM_Face.h

@@ -29,6 +29,7 @@ class FM_Face {
     }
   inline void AddEdge(FM_Edge* edge) { _edge.push_back(edge); }
   inline int GetNumEdges() { return _edge.size(); }
+  inline Patch* GetPatch() { return _patch; }
   inline FM_Edge* GetEdge(int i) { return _edge[i]; }
 
   void F(double u, double v, double &x, double &y, double &z);

contrib/FourierModel/FM_Info.cpp

@@ -4,6 +4,8 @@ PatchInfo::PatchInfo() : tag(-1)
 {
   uMin = vMin = 0.;
   uMax = vMax = 1.;
+  derivative = 3;
+  recompute = 1;
 }
 
 OverlapInfo::OverlapInfo()

contrib/FourierModel/FM_Info.h

@@ -6,15 +6,20 @@
 
 class PatchInfo {
  public:
-  char type[16];
   int tag;
-  int projectionSurfaceTag;
+  int derivative;
+  int recompute;
   int nM[2];
   int nModes[2];
-  int periodic[2];
   double uMin, uMax, vMin, vMax;
   int hardEdge[4];
-  std::vector<std::vector<std::complex<double> > > coeff;
+  std::vector<std::vector<std::complex<double> > > coeffFourier;
+  std::vector<std::vector<std::complex<double> > > coeffCheby;
+  std::vector<std::vector<std::complex<double> > > coeffDerivU;
+  std::vector<std::vector<std::complex<double> > > coeffDerivV;
+  std::vector<std::vector<std::complex<double> > > coeffDerivUU;
+  std::vector<std::vector<std::complex<double> > > coeffDerivUV;
+  std::vector<std::vector<std::complex<double> > > coeffDerivVV;
 
   PatchInfo();
   virtual ~PatchInfo() {}

contrib/FourierModel/FM_Reader.cpp

@@ -6,156 +6,145 @@
 FM_Reader::FM_Reader(const char* fn)
 {
   char c;
-  char continuation[16] = "continuation";
+  char Exact[8] = "Exact";
   char cylinder[16] = "cylinder";
   char revolvedParabola[32] = "revolvedParabola";
+
   std::ifstream InputFile(fn);
   if (!InputFile) {
     Msg::Info("Failed to open input file.");
     exit(EXIT_FAILURE);
   }
-  InputFile >> _nPS;
-  for (int i=0;i<_nPS;i++) {
-    char projection[32];
-    InputFile >> projection;
+  //InputFile >> _nPatches;
+  //for (unsigned int i = 0; i < _nPatches; i++) {
+  unsigned int i = 0;
+  while (!InputFile.eof()) {
+    std::cout << "here\n";
+    char psName[32];
+    InputFile >> psName;
+    std::cout << psName;
+    int psTag;
+    InputFile >> psTag;
     double origin[3];
     InputFile >> origin[0] >> origin[1] >> origin[2];
-    double normal[3];
-    InputFile >> normal[0] >> normal[1] >> normal[2];
-    double tangent[3];
-    InputFile >> tangent[0] >> tangent[1] >> tangent[2];
+    double E0[3];
+    InputFile >> E0[0] >> E0[1] >> E0[2];
+    double E1[3];
+    InputFile >> E1[0] >> E1[1] >> E1[2];
     double scale[3];
     InputFile >> scale[0] >> scale[1] >> scale[2];
-    if (!strcmp(projection,cylinder))
+    int psNumParams;
+    InputFile >> psNumParams;
+    std::vector<double> psParams;
+    for (unsigned int j = 0; j < psNumParams; j++) {
+      double tmp;
+      InputFile >> tmp;
+      psParams.push_back(tmp);
+    }
+    if (!strcmp(psName,cylinder))
       _ps.push_back
-	(new CylindricalProjectionSurface(i,origin,normal,tangent,scale));
-    else if (!strcmp(projection,revolvedParabola)) {
-      double R, K[2];
-      InputFile >> R;
-      InputFile >> K[0] >> K[1];
+	(new CylindricalProjectionSurface
+	 (psTag,origin,E0,E1,scale,psParams[0],psParams[1]));
+    else if (!strcmp(psName,revolvedParabola)) {
+      double K[2];
+      K[0] = psParams[1]; K[1] = psParams[2];
       _ps.push_back(new RevolvedParabolaProjectionSurface
-		    (i,origin,normal,tangent,scale,R,K));
+		    (i,origin,E0,E1,scale,psParams[0],K));
     }
     else {
       _ps.push_back
-	(new CylindricalProjectionSurface(i,origin,normal,tangent,scale));
+	(new CylindricalProjectionSurface(i,origin,E0,E1,scale));
       Msg::Error("Unknown projection surface. Replaced by Cylinder...");
     }
-  }
-
-  InputFile >> _nPatches;
-
-  for (int i=0;i<_nPatches;i++) {
+    InputFile >> psName;
     _patchList.push_back(new PatchInfo);
-    InputFile >> _patchList[i]->type;
-    InputFile >> _patchList[i]->projectionSurfaceTag;
-    if (!strcmp(_patchList[i]->type,continuation)) {
-      InputFile >> _patchList[i]->periodic[0] >> _patchList[i]->periodic[1];
-      InputFile >> _patchList[i]->uMin >> _patchList[i]->uMax;
-      InputFile >> _patchList[i]->vMin >> _patchList[i]->vMax;
+    InputFile >> _patchList[i]->tag;
+    InputFile >> _patchList[i]->tag;
+    InputFile >> _patchList[i]->uMin >> _patchList[i]->uMax;
+    InputFile >> _patchList[i]->vMin >> _patchList[i]->vMax;
+    if (strcmp(psName,Exact)) {
       InputFile >> _patchList[i]->hardEdge[0] >> _patchList[i]->hardEdge[1] >>
 	_patchList[i]->hardEdge[2] >> _patchList[i]->hardEdge[3];
-      InputFile >> _patchList[i]->nM[0] >> _patchList[i]->nM[1];
       InputFile >> _patchList[i]->nModes[0] >> _patchList[i]->nModes[1];
-      _patchList[i]->coeff.resize(_patchList[i]->nModes[0]);
+      _patchList[i]->coeffFourier.resize(_patchList[i]->nModes[0]);
       for (int j=0;j<_patchList[i]->nModes[0];j++) {
-	_patchList[i]->coeff[j].resize(_patchList[i]->nModes[1]);
+	_patchList[i]->coeffFourier[j].resize(_patchList[i]->nModes[1]);
 	for (int k=0;k<_patchList[i]->nModes[1];k++) {
 	  double realCoeff, imagCoeff;
 	  InputFile >> realCoeff >> imagCoeff;
-	  _patchList[i]->coeff[j][k] = 
+	  _patchList[i]->coeffFourier[j][k] = 
 	    std::complex<double>(realCoeff,imagCoeff);
 	}
       }
+      InputFile >> _patchList[i]->nM[0] >> _patchList[i]->nM[1];
+      InputFile >> _patchList[i]->recompute;
+      if ((_patchList[i]->derivative) && (!_patchList[i]->recompute)) {
+	_patchList[i]->coeffCheby.resize(_patchList[i]->nM[0]);
+	for (int j=0;j<_patchList[i]->nM[0];j++) {
+	  _patchList[i]->coeffCheby[j].resize(_patchList[i]->nM[1]);
+	  for (int k=0;k<_patchList[i]->nM[1];k++) {
+	    double realCoeff, imagCoeff;
+	    InputFile >> realCoeff >> imagCoeff;
+	    _patchList[i]->coeffCheby[j][k] = 
+	      std::complex<double>(realCoeff,imagCoeff);
+	  }
+	}
+	_patchList[i]->coeffDerivU.resize(_patchList[i]->nM[0]);
+	for (int j=0;j<_patchList[i]->nM[0];j++) {
+	  _patchList[i]->coeffDerivU[j].resize(_patchList[i]->nM[1]);
+	  for (int k=0;k<_patchList[i]->nM[1];k++) {
+	    double realCoeff, imagCoeff;
+	    InputFile >> realCoeff >> imagCoeff;
+	    _patchList[i]->coeffDerivU[j][k] = 
+	      std::complex<double>(realCoeff,imagCoeff);
+	  }
+	}
+	_patchList[i]->coeffDerivV.resize(_patchList[i]->nM[0]);
+	for (int j=0;j<_patchList[i]->nM[0];j++) {
+	  _patchList[i]->coeffDerivV[j].resize(_patchList[i]->nM[1]);
+	  for (int k=0;k<_patchList[i]->nM[1];k++) {
+	    double realCoeff, imagCoeff;
+	    InputFile >> realCoeff >> imagCoeff;
+	    _patchList[i]->coeffDerivV[j][k] = 
+	      std::complex<double>(realCoeff,imagCoeff);
+	  }
+	}
+	_patchList[i]->coeffDerivUU.resize(_patchList[i]->nM[0]);
+	for (int j=0;j<_patchList[i]->nM[0];j++) {
+	  _patchList[i]->coeffDerivUU[j].resize(_patchList[i]->nM[1]);
+	  for (int k=0;k<_patchList[i]->nM[1];k++) {
+	    double realCoeff, imagCoeff;
+	    InputFile >> realCoeff >> imagCoeff;
+	    _patchList[i]->coeffDerivUU[j][k] = 
+	      std::complex<double>(realCoeff,imagCoeff);
+	  }
+	}
+	_patchList[i]->coeffDerivUV.resize(_patchList[i]->nM[0]);
+	for (int j=0;j<_patchList[i]->nM[0];j++) {
+	  _patchList[i]->coeffDerivUV[j].resize(_patchList[i]->nM[1]);
+	  for (int k=0;k<_patchList[i]->nM[1];k++) {
+	    double realCoeff, imagCoeff;
+	    InputFile >> realCoeff >> imagCoeff;
+	    _patchList[i]->coeffDerivUV[j][k] = 
+	      std::complex<double>(realCoeff,imagCoeff);
+	  }
+	}
+	_patchList[i]->coeffDerivVV.resize(_patchList[i]->nM[0]);
+	for (int j=0;j<_patchList[i]->nM[0];j++) {
+	  _patchList[i]->coeffDerivVV[j].resize(_patchList[i]->nM[1]);
+	  for (int k=0;k<_patchList[i]->nM[1];k++) {
+	    double realCoeff, imagCoeff;
+	    InputFile >> realCoeff >> imagCoeff;
+	    _patchList[i]->coeffDerivVV[j][k] = 
+	      std::complex<double>(realCoeff,imagCoeff);
+	  }
+	}
+      }
     }
+    _patch.push_back(new ContinuationPatch(_patchList[i], _ps[i]));
+    i++;
   }
-
-  InputFile >> _nIntersections;
-
-  for (int i=0;i<_nIntersections;i++) {
-    _intersectionList.push_back(new IntersectionInfo);
-    InputFile >> _intersectionList[i]->tag;
-    InputFile >> _intersectionList[i]->SP[0] >> _intersectionList[i]->SP[1] 
-	      >> _intersectionList[i]->SP[2];
-    InputFile >> _intersectionList[i]->EP[0] >> _intersectionList[i]->EP[1] 
-	      >> _intersectionList[i]->EP[2];
-    InputFile >> _intersectionList[i]->intersectingPatches[0].patchTag;
-    InputFile >> _intersectionList[i]->intersectingPatches[1].patchTag;
-    InputFile >> _intersectionList[i]->along;
-    if ((_intersectionList[i]->intersectingPatches[0].patchTag < 0) ||
-	(_intersectionList[i]->intersectingPatches[1].patchTag < 0))
-      InputFile >> _intersectionList[i]->edgeInfo.edgeType >>
-	_intersectionList[i]->edgeInfo.constValue >>
-	_intersectionList[i]->edgeInfo.startValue >>
-	_intersectionList[i]->edgeInfo.endValue >>
-	_intersectionList[i]->edgeInfo.acrossDiscontinuity;
-  }
-  _overlapChart.resize(_nPatches);
-  for (int i=0;i<_nPatches;i++) {
-    _overlapChart[i].resize(_nPatches);
-    for (int j=0;j<_nPatches;j++) {
-      _overlapChart[i][j] = new OverlapInfo;
-      InputFile >> _overlapChart[i][j]->doesIntersect >> 
-	_overlapChart[i][j]->xMin >> _overlapChart[i][j]->xMax >>
-	_overlapChart[i][j]->yMin >> _overlapChart[i][j]->yMax >>
-	_overlapChart[i][j]->zMin >> _overlapChart[i][j]->zMax >>
-	_overlapChart[i][j]->psTag;
-    }
-  }
-
-  _curve.resize(_nIntersections,0);
-
-  for (int i=0;i<_nPatches;i++) {
-    _patchList[i]->tag = i;
-    PatchInfo* PI = _patchList[i];
-
-    if (!strcmp(PI->type,continuation))
-      _patch.push_back
-	(new ContinuationPatch(PI,_ps[PI->projectionSurfaceTag],3));
-    else
-      _patch.push_back(new ExactPatch(PI,_ps[PI->projectionSurfaceTag]));
-  }
-
-  _blendOperator = new BlendOperator(_patch,_ps,_overlapChart);
-
-  for (int i=0;i<_nPatches;i++)
-    _blendedPatch.push_back(new BlendedPatch(_patch[i],_blendOperator));
-
-  for (int i=0;i<_nIntersections;i++) {
-    IntersectionInfo* II = _intersectionList[i];
-    _curve[II->tag] = new IntersectionCurve(II,_patch);
-  }
-
-  InputFile >> _nVertices;
-  for (int i=0;i<_nVertices;i++) {
-    double x,y,z;
-    InputFile >> x >> y >> z;
-    _vertex.push_back(new FM_Vertex(i,x,y,z));
-  }
-
-  InputFile >> _nEdges;
-  for (int i=0;i<_nEdges;i++) {
-    int edgeTag, svTag, evTag;
-    InputFile >> edgeTag >> svTag >> evTag;
-    if (edgeTag < 0)
-      _edge.push_back(new FM_Edge(i,0,_vertex[svTag],_vertex[evTag]));
-    else
-      _edge.push_back(new FM_Edge(i,GetCurve(edgeTag),
-				  _vertex[svTag],_vertex[evTag]));
-  }
-
-  InputFile >> _nFaces;
-  for (int i=0;i<_nFaces;i++) {
-    int faceTag, nEdges;
-    InputFile >> faceTag;
-    _face.push_back(new FM_Face(i,GetPatch(faceTag)));
-    InputFile >> nEdges;
-    for (int j=0;j<nEdges;j++) {
-      int edgeTag;
-      InputFile >> edgeTag;
-      _face[i]->AddEdge(_edge[edgeTag]);
-    }
-  }
+  _nPatches = _patch.size();
 }
 
 Patch* FM_Reader::GetPatch(int tag)
@@ -165,18 +154,6 @@ Patch* FM_Reader::GetPatch(int tag)
       return _patch[i];
 }
 
-Curve* FM_Reader::GetCurve(int tag)
-{
-  Curve* curve = 0;
-  for (int i=0;i<_curve.size();i++) {
-    if (_curve[i]->GetTag() == tag) {
-      curve = _curve[i];
-      break;
-    }
-  }
-  return curve;
-}
-
 ProjectionSurface* FM_Reader::GetProjectionSurface(int tag)
 {
   ProjectionSurface* ps = 0;
@@ -188,10 +165,3 @@ ProjectionSurface* FM_Reader::GetProjectionSurface(int tag)
   }
   return ps;
 }
-
-BlendedPatch* FM_Reader::GetBlendedPatch(int tag)
-{
-  for (int i=0;i<_blendedPatch.size();i++)
-    if (_blendedPatch[i]->GetTag() == tag)
-      return _blendedPatch[i];
-}

contrib/FourierModel/FM_Reader.h

@@ -17,76 +17,25 @@
 
 class FM_Reader {
  private:
-  int _nVertices;
-  int _nEdges;
-  int _nFaces;
   int _nPatches;
-  int _nPS;
-  int _nIntersections;
   std::vector<PatchInfo*> _patchList;
-  std::vector<IntersectionInfo*> _intersectionList;
-  std::vector<std::vector<OverlapInfo*> > _overlapChart;
-  std::vector<FM_Vertex*> _vertex;
-  std::vector<FM_Edge*> _edge;
-  std::vector<FM_Face*> _face;
   std::vector<Patch*> _patch;
-  std::vector<Curve*> _curve;
   std::vector<ProjectionSurface*> _ps;
 
-  BlendOperator* _blendOperator;
-  std::vector<BlendedPatch*> _blendedPatch;
-
  public:
   FM_Reader(const char* fn);
 
   virtual ~FM_Reader() {}
 
-  inline int GetNumVertices
-    () { return _nVertices; }
-
-  inline int GetNumEdges
-    () { return _nEdges; }
-
-  inline int GetNumFaces
-    () { return _nFaces; }
-
   inline int GetNumPatches
     () { return _nPatches; }
 
-  inline int GetNumIntersections
-    () { return _nIntersections; }
-
   inline std::vector<PatchInfo*> GetPatchList
     () { return _patchList; }
 
-  inline std::vector<IntersectionInfo*> GetIntersectionList
-    () { return _intersectionList; }
-
-  inline FM_Vertex* GetVertex
-    (int i) { return _vertex[i]; }
-
-  inline FM_Edge* GetEdge
-    (int i) { return _edge[i]; }
-
-  inline FM_Face* GetFace
-    (int i) { return _face[i]; }
-
-  inline bool GetOverlapInfo
-    (int i, int j) 
-    { if (_overlapChart[i][j]->doesIntersect) return true; else return false; }
-
-  inline BlendOperator* GetBlendOperator
-    () { return _blendOperator; }
-
   Patch* GetPatch
     (int tag);
 
-  BlendedPatch* GetBlendedPatch
-    (int tag);
-
-  Curve* GetCurve
-    (int tag);
-
   ProjectionSurface* GetProjectionSurface
     (int tag);
 };

contrib/FourierModel/FPatch.cpp

@@ -855,3 +855,78 @@ double  FPatch::GetPou(double u, double v)
 
   return pouU * pouV;
 }
+
+void FPatch::Export(FILE *fp)
+{
+  double x,y,z;
+
+  fprintf(fp, "\n");
+  fprintf(fp, "%s\n", _ps->GetName().c_str());
+  fprintf(fp, "%d\n", _ps->GetTag());
+  _ps->GetOrigin(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  _ps->GetE0(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  _ps->GetE1(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  _ps->GetScale(x,y,z);
+  fprintf(fp, "%g %g %g\n", x, y, z);
+  fprintf(fp, "%d\n", _ps->GetNumParameters());
+  for (int i = 0; i < _ps->GetNumParameters(); i++)
+    fprintf(fp, "%g\n",_ps->GetParameter(i));
+
+  fprintf(fp, "%s\n", "FPatch");
+  fprintf(fp, "%d\n", _tag);
+  fprintf(fp, "%d\n", _derivative);
+  fprintf(fp, "%g %g\n", _uMin, _uMax);
+  fprintf(fp, "%g %g\n", _vMin, _vMax);
+  fprintf(fp, "%d %d %d %d\n", _hardEdge[0], _hardEdge[1], _hardEdge[2],
+	  _hardEdge[3]);
+  fprintf(fp, "%d %d\n", _uModes, _vModes);
+  for (int j = 0; j < _uModes; j++) {
+    for (int k = 0; k < _vModes; k++) {
+      fprintf(fp, "%g %g\n", _coeffOriginalData[j][k].real(),
+	      _coeffOriginalData[j][k].imag());
+    }
+  }
+  fprintf(fp, "%d %d\n", _uM, _vM);
+  fprintf(fp, "%d\n", 0);
+  if (_derivative) {
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffData[j][k].real(), 
+		_coeffData[j][k].imag());
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffDerivU[j][k].real(), 
+		_coeffDerivU[j][k].imag());
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffDerivV[j][k].real(), 
+		_coeffDerivV[j][k].imag());
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffDerivUU[j][k].real(), 
+		_coeffDerivUU[j][k].imag());
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffDerivUV[j][k].real(), 
+		_coeffDerivUV[j][k].imag());
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	fprintf(fp, "%g %g\n", _coeffDerivVV[j][k].real(), 
+		_coeffDerivVV[j][k].imag());
+      }
+    }
+  }
+}

contrib/FourierModel/FPatch.h

@@ -78,6 +78,7 @@ class FPatch : public Patch {
 
   // Abstract functions of Patch
 
+  virtual void Export(FILE *fp);
   virtual double GetPou(double u, double v);
   virtual void F(double u, double v, double &x, double &y, double &z);
   virtual bool Inverse(double x,double y,double z,double &u,double &v);

contrib/FourierModel/Makefile

@@ -5,6 +5,8 @@ LIB = ../../lib/libGmshFourierModel.a
 CFLAGS = ${OPTIM} ${FLAGS}
 
 SRC = ProjectionSurface.cpp \
+	PlaneProjectionSurface.cpp \
+	ParaboloidProjectionSurface.cpp \
 	CylindricalProjectionSurface.cpp \
 	RevolvedParabolaProjectionSurface.cpp \
       Patch.cpp \
@@ -49,6 +51,9 @@ depend:
 
 # DO NOT DELETE THIS LINE
 ProjectionSurface.o: ProjectionSurface.cpp ProjectionSurface.h
+PlaneProjectionSurface.o: PlaneProjectionSurface.cpp PlaneProjectionSurface.h
+ParaboloidProjectionSurface.o: ParaboloidProjectionSurface.cpp \
+  ParaboloidProjectionSurface.h
 CylindricalProjectionSurface.o: CylindricalProjectionSurface.cpp \
   CylindricalProjectionSurface.h ProjectionSurface.h
 RevolvedParabolaProjectionSurface.o:  \

contrib/FourierModel/Patch.cpp

@@ -2,7 +2,7 @@
 #include "Patch.h"
 
 Patch::Patch() :_ps(0), _uMin(0.), _uMax(1.), _vMin(0.), _vMax(1.),
-		_periodicityU(0), _periodicityV(0), _derivative(0), 
+		_periodicityU(0), _periodicityV(0), _derivative(3), 
 		_tag(-1) {}
 
 void Patch::GetNormal(double u, double v, double &x, double &y, double &z)