From 87117cc064d6a2084b30c12a1e4fdddb41d59321 Mon Sep 17 00:00:00 2001
From: Akash Anand <akasha@iitk.ac.in>
Date: Sat, 12 May 2007 22:45:14 +0000
Subject: [PATCH] *** empty log message ***
---
Geo/FEdge.cpp | 68 +--
Geo/FEdge.h | 3 -
Geo/FFace.cpp | 2 +-
Mesh/meshGEdge.cpp | 3 +-
contrib/FourierModel/ContinuationPatch.cpp | 4 +-
contrib/FourierModel/Curve.cpp | 463 ++++++++++++---------
contrib/FourierModel/Curve.h | 1 +
contrib/FourierModel/FM_Edge.cpp | 24 +-
contrib/FourierModel/FM_Info.h | 3 +-
contrib/FourierModel/FM_Reader.cpp | 17 +-
10 files changed, 310 insertions(+), 278 deletions(-)
diff --git a/Geo/FEdge.cpp b/Geo/FEdge.cpp
index 34ed1e4667..82b9f7eabb 100644
--- a/Geo/FEdge.cpp
+++ b/Geo/FEdge.cpp
@@ -4,19 +4,13 @@
#if defined(HAVE_FOURIER_MODEL)
FEdge::FEdge(GModel *model, FM_Edge* edge_, int tag, GVertex *v0, GVertex *v1)
- : GEdge(model, tag, v0, v1), edge(edge_), face(0), edgeNum(-1)
+ : GEdge(model, tag, v0, v1), edge(edge_)
{
//meshAttributes.Method = TRANSFINI;
//meshAttributes.coeffTransfinite = 1.;
//meshAttributes.nbPointsTransfinite = 10;
}
-FEdge::FEdge(GModel *model, FM_Face* face_, int edgeNum_, int tag, GVertex *v0,
- GVertex *v1)
- : GEdge(model, tag, v0, v1), edge(0), face(face_), edgeNum(edgeNum_)
-{
-}
-
Range<double> FEdge::parBounds(int i) const
{
return(Range<double>(0.,1.));
@@ -25,73 +19,23 @@ Range<double> FEdge::parBounds(int i) const
GPoint FEdge::point(double p) const
{
double x, y, z;
- if (edge)
- edge->F(p,x,y,z);
- else {
- if (edgeNum == 0) {
- double p0, p1, q0, q1;
- face->Inverse(v0->x(),v0->y(),v0->z(),p0,q0);
- face->Inverse(v1->x(),v1->y(),v1->z(),p1,q1);
- face->F(p0+p*(p1-p0),0,x,y,z);
- }
- else if (edgeNum == 1) {
- double p0, p1, q0, q1;
- face->Inverse(v0->x(),v0->y(),v0->z(),q0,p0);
- face->Inverse(v1->x(),v1->y(),v1->z(),q1,p1);
- face->F(1.,p0+p*(p1-p0),x,y,z);
- }
- else if (edgeNum == 2) {
- double p0, p1, q0, q1;
- face->Inverse(v0->x(),v0->y(),v0->z(),p0,q0);
- face->Inverse(v1->x(),v1->y(),v1->z(),p1,q1);
- face->F(p0+p*(p1-p0),1.,x,y,z);
- }
- else if (edgeNum == 3) {
- double p0, p1, q0, q1;
- face->Inverse(v0->x(),v0->y(),v0->z(),q0,p0);
- face->Inverse(v1->x(),v1->y(),v1->z(),q1,p1);
- face->F(0.,p0+p*(p1-p0),x,y,z);
- }
- else
- Msg(INFO,"Invalid edge number.");
- }
+ edge->F(p,x,y,z);
+
return GPoint(x,y,z);
}
double FEdge::parFromPoint(const SPoint3 &pt) const
{
double p;
- if (edge)
- edge->Inverse(pt.x(),pt.y(),pt.z(),p);
- else {
- if ((edgeNum == 0) || (edgeNum == 2)) {
- double p0, p1, q0, q1, q;
- face->Inverse(v0->x(),v0->y(),v0->z(),p0,q0);
- face->Inverse(v1->x(),v1->y(),v1->z(),p1,q1);
- face->Inverse(pt.x(),pt.y(),pt.z(),p,q);
- p = p0 + p * (p1 - p0);
- }
- else if ((edgeNum == 1) || (edgeNum == 3)) {
- double p0, p1, q0, q1, q;
- face->Inverse(v0->x(),v0->y(),v0->z(),q0,p0);
- face->Inverse(v1->x(),v1->y(),v1->z(),q1,p1);
- face->Inverse(pt.x(),pt.y(),pt.z(),q,p);
- p = p0 + p * (p1 - p0);
- }
- else
- Msg(INFO,"Invalid edge number.");
- }
+ edge->Inverse(pt.x(),pt.y(),pt.z(),p);
+
return p;
}
SVector3 FEdge::firstDer(double par) const
{
double x,y,z;
- if (edge)
- edge->Dfdt(par,x,y,z);
- else {
- x = y = z = 0.;
- }
+ edge->Dfdt(par,x,y,z);
return SVector3(x,y,z);
}
diff --git a/Geo/FEdge.h b/Geo/FEdge.h
index 6ff0bb9987..b11b79d4d5 100644
--- a/Geo/FEdge.h
+++ b/Geo/FEdge.h
@@ -15,12 +15,9 @@
class FEdge : public GEdge {
protected:
FM_Edge* edge;
- FM_Face* face;
int edgeNum;
public:
FEdge(GModel *model, FM_Edge* edge_, int tag, GVertex *v0, GVertex *v1);
- FEdge(GModel *model, FM_Face* face_, int edgeNum_, int tag, GVertex *v0,
- GVertex *v1);
virtual ~FEdge() {}
double period() const { throw ; }
virtual bool periodic(int dim=0) const { return false; }
diff --git a/Geo/FFace.cpp b/Geo/FFace.cpp
index 5040ad75c0..043191f86d 100644
--- a/Geo/FFace.cpp
+++ b/Geo/FFace.cpp
@@ -58,7 +58,7 @@ int FFace::containsParam(const SPoint2 &pt) const
if(pt[1] < 0. - tol || pt[1] > 1. + tol) return 0;
return 1;
}
-
+
SVector3 FFace::normal(const SPoint2 ¶m) const
{
double x,y,z;
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 61824ec7c4..81c23aa335 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.38 2007-05-02 07:59:27 geuzaine Exp $
+// $Id: meshGEdge.cpp,v 1.39 2007-05-12 22:45:14 anand Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -227,6 +227,7 @@ void meshGEdge::operator() (GEdge *ge)
a = Integration(ge, t_begin, t_end, F_Lc, Points, 1.e-8);
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
}
+
const double b = a / (double)(N - 1);
// if the curve is periodic and if the begin vertex is identical to the end vertex
diff --git a/contrib/FourierModel/ContinuationPatch.cpp b/contrib/FourierModel/ContinuationPatch.cpp
index 8018987703..7f11b0e988 100644
--- a/contrib/FourierModel/ContinuationPatch.cpp
+++ b/contrib/FourierModel/ContinuationPatch.cpp
@@ -166,8 +166,8 @@ std::complex<double> ContinuationPatch::
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.);
- }
+ "not in [%g,%g]x[%g,%g]", u, v, 0., 1., 0., 1.);
+ }
// Interpolate to find value at (u,v)
for(int j = 0; j < _uModes; j++){
diff --git a/contrib/FourierModel/Curve.cpp b/contrib/FourierModel/Curve.cpp
index 81191ac27f..fbe26d710a 100644
--- a/contrib/FourierModel/Curve.cpp
+++ b/contrib/FourierModel/Curve.cpp
@@ -2,228 +2,307 @@
#include "Message.h"
#include "Curve.h"
-Curve::Curve(IntersectionInfo* II, std::vector<Patch*> patches) : _II(II)
+Curve::Curve(IntersectionInfo* II, std::vector<Patch*> patches)
+ : _II(II)
{
_h = 1.e-6;
_tol = 1.e-12;
- for (int i=0;i<patches.size();i++) {
- if (patches[i]->GetTag() == _II->intersectingPatches[0].patchTag) {
- _patch0 = patches[i];
- break;
+ if (_II->intersectingPatches[0].patchTag < 0)
+ _patch0 = 0;
+ else {
+ for (int i=0;i<patches.size();i++) {
+ if (patches[i]->GetTag() == _II->intersectingPatches[0].patchTag) {
+ _patch0 = patches[i];
+ break;
+ }
}
}
- for (int i=0;i<patches.size();i++)
- if (patches[i]->GetTag() == _II->intersectingPatches[1].patchTag) {
- _patch1 = patches[i];
- break;
- }
- /*
- int n1=64;
- int n2=64;
- double h1 = 1000./(n1-1);
- double h2 = 1000./(n2-1);
- for (int j=0;j<n1;j++) {
- for (int k=0;k<n2;k++) {
- double u = j*h1;
- double v = k*h2;
- double x, y, z;
- _patch1->F(u,v,x,y,z);
- printf("x(%d,%d) = %g; y(%d,%d) = %g; z(%d,%d) = %g;\n",
- j+1,k+1,x,j+1,k+1,y,j+1,k+1,z);
+ if (_II->intersectingPatches[1].patchTag)
+ _patch1 = 0;
+ else {
+ for (int i=0;i<patches.size();i++) {
+ if (patches[i]->GetTag() == _II->intersectingPatches[1].patchTag) {
+ _patch1 = patches[i];
+ break;
+ }
}
}
- */
- double u,v,x,y,z;
- x = _II->SP[0]; y = _II->SP[1]; z = _II->SP[2];
- _patch0->Inverse(x,y,z,u,v);
- _SP0[0] = u;
- _SP0[1] = v;
- _patch1->Inverse(x,y,z,u,v);
- _SP1[0] = u;
- _SP1[1] = v;
- //printf("SPx = %g; SPy = %g; SPz = %g;\n",x,y,z);
-
- x = _II->EP[0]; y = _II->EP[1]; z = _II->EP[2];
- _patch0->Inverse(x,y,z,u,v);
- _EP0[0] = u;
- _EP0[1] = v;
- _patch1->Inverse(x,y,z,u,v);
- _EP1[0] = u;
- _EP1[1] = v;
- //printf("EPx = %g; EPy = %g; EPz = %g;\n",x,y,z);
+ if ((_patch0) && (_patch1)) {
+ double u,v,x,y,z;
+ x = _II->SP[0]; y = _II->SP[1]; z = _II->SP[2];
+ _patch0->Inverse(x,y,z,u,v);
+ _SP0[0] = u;
+ _SP0[1] = v;
+ _patch1->Inverse(x,y,z,u,v);
+ _SP1[0] = u;
+ _SP1[1] = v;
+ //printf("SPx = %g; SPy = %g; SPz = %g;\n",x,y,z);
+
+ x = _II->EP[0]; y = _II->EP[1]; z = _II->EP[2];
+ _patch0->Inverse(x,y,z,u,v);
+ _EP0[0] = u;
+ _EP0[1] = v;
+ _patch1->Inverse(x,y,z,u,v);
+ _EP1[0] = u;
+ _EP1[1] = v;
+ //printf("EPx = %g; EPy = %g; EPz = %g;\n",x,y,z);
+ }
+ else if (_patch0) {
+ double u,v,x,y,z;
+ x = _II->SP[0]; y = _II->SP[1]; z = _II->SP[2];
+ _patch0->Inverse(x,y,z,u,v);
+ _SP0[0] = u;
+ _SP0[1] = v;
+ x = _II->EP[0]; y = _II->EP[1]; z = _II->EP[2];
+ _patch0->Inverse(x,y,z,u,v);
+ _EP0[0] = u;
+ _EP0[1] = v;
+ }
+ else if (_patch1) {
+ double u,v,x,y,z;
+ x = _II->SP[0]; y = _II->SP[1]; z = _II->SP[2];
+ _patch1->Inverse(x,y,z,u,v);
+ _SP1[0] = u;
+ _SP1[1] = v;
+ x = _II->EP[0]; y = _II->EP[1]; z = _II->EP[2];
+ _patch1->Inverse(x,y,z,u,v);
+ _EP1[0] = u;
+ _EP1[1] = v;
+ }
+ else
+ Msg::Error("No patches to find the intersection curve.");
}
void Curve::F(double t, double &x, double &y, double &z)
{
- double u0, v0, u1, v1;
- //Msg::Info("patch0 : %d",_patch0->GetTag());
- //Msg::Info("patch1 : %d",_patch1->GetTag());
- //Msg::Info("%g,%g,%g",_II->SP[0],_II->SP[1],_II->SP[2]);
- //Msg::Info("%g,%g,%g",_II->EP[0],_II->EP[1],_II->EP[2]);
- //Msg::Info("%g,%g,%g,%g",_SP0[0],_SP0[1],_EP0[0],_EP0[1]);
- //Msg::Info("%g,%g,%g,%g",_SP1[0],_SP1[1],_EP1[0],_EP1[1]);
-
- if ((std::abs(_SP0[0]-_EP0[0])<1.e-12) &&
- (_patch0->_PI->periodic[0] == 1)) {
- u0 = _SP0[0] + t * (1. + _EP0[0] - _SP0[0]);
- if (u0 > 1.)
+ if ((_patch0) && (_patch1)) {
+ double u0, v0, u1, v1;
+ //printf("patch0 : %d\n",_patch0->GetTag());
+ //printf("patch1 : %d\n",_patch1->GetTag());
+ //printf("%g,%g,%g\n",_II->SP[0],_II->SP[1],_II->SP[2]);
+ //printf("%g,%g,%g\n",_II->EP[0],_II->EP[1],_II->EP[2]);
+ //printf("%g,%g,%g,%g\n",_SP0[0],_SP0[1],_EP0[0],_EP0[1]);
+ //printf("%g,%g,%g,%g\n",_SP1[0],_SP1[1],_EP1[0],_EP1[1]);
+
+ if (_patch0->_PI->periodic[0] == 1) {
+ if (std::abs(_SP0[0]-_EP0[0])<1.e-12)
+ u0 = _SP0[0] + t * (1. + _EP0[0] - _SP0[0]);
+ else
+ u0 = _SP0[0] + t * (_EP0[0] - _SP0[0]);
u0 -= floor(u0);
- }
- else
- u0 = _SP0[0] + t * (_EP0[0] - _SP0[0]);
-
- if ((std::abs(_SP0[1]-_EP0[1])<1.e-12) &&
- (_patch0->_PI->periodic[1] == 1)) {
- v0 = _SP0[1] + t * (1. + _EP0[1] - _SP0[1]);
- if (u0 > 1.)
+ }
+ else
+ u0 = _SP0[0] + t * (_EP0[0] - _SP0[0]);
+
+ if (_patch0->_PI->periodic[1] == 1) {
+ if (std::abs(_SP0[1]-_EP0[1])<1.e-12)
+ v0 = _SP0[1] + t * (1. + _EP0[1] - _SP0[1]);
+ else
+ v0 = _SP0[1] + t * (_EP0[1] - _SP0[1]);
v0 -= floor(v0);
- }
- else
- v0 = _SP0[1] + t * (_EP0[1] - _SP0[1]);
-
- if ((std::abs(_SP1[0]-_EP1[0])<1.e-12) &&
- (_patch1->_PI->periodic[0] == 1)) {
- u1 = _SP1[0] + t * (1. + _EP1[0] - _SP1[0]);
- if (u1 > 1.)
+ }
+ else
+ v0 = _SP0[1] + t * (_EP0[1] - _SP0[1]);
+
+ if (_patch1->_PI->periodic[0] == 1) {
+ if (std::abs(_SP1[0]-_EP1[0])<1.e-12)
+ u1 = _SP1[0] + t * (1. + _EP1[0] - _SP1[0]);
+ else
+ u1 = _SP1[0] + t * (_EP1[0] - _SP1[0]);
u1 -= floor(u1);
- }
- else {
- u1 = _SP1[0] + t * (_EP1[0] - _SP1[0]);
- }
-
- if ((std::abs(_SP1[1]-_EP1[1])<1.e-12) &&
- (_patch1->_PI->periodic[1] == 1)) {
- v1 = _SP1[1] + t * (1. + _EP1[1] - _SP1[1]);
- if (u1 > 1.)
+ }
+ else
+ u1 = _SP1[0] + t * (_EP1[0] - _SP1[0]);
+
+ if (_patch1->_PI->periodic[1] == 1) {
+ if (std::abs(_SP1[1]-_EP1[1])<1.e-12)
+ v1 = _SP1[1] + t * (1. + _EP1[1] - _SP1[1]);
+ else
+ v1 = _SP1[1] + t * (_EP1[1] - _SP1[1]);
v1 -= floor(v1);
- }
- else
- v1 = _SP1[1] + t * (_EP1[1] - _SP1[1]);
-
- //Msg::Info("%g,%g,%g,%g",u0,v0,u1,v1);
-
- double x0, y0, z0;
- _patch0->F(u0,v0,x0,y0,z0);
-
- //Msg::Info("%g,%g,%g",x0,y0,z0);
-
- double x1, y1, z1;
- _patch1->F(u1,v1,x1,y1,z1);
-
- //Msg::Info("%g,%g,%g",x1,y1,z1);
-
- double r,u,v,rPlus,uPlus,vPlus;
-
- r = u0;
- u = u1;
- v = v1;
-
- double F[3];
- F[0] = x1 - x0;
- F[1] = y1 - y0;
- F[2] = z1 - z0;
-
- while ((std::abs(F[0])>_tol) || (std::abs(F[1])>_tol) ||
- (std::abs(F[2])>_tol)) {
- rPlus = r + _h;
- uPlus = u + _h;
- vPlus = v + _h;
- if (_patch0->_PI->periodic[0] == 1)
- rPlus -= std::floor(rPlus);
- if (_patch1->_PI->periodic[0] == 1)
- uPlus -= std::floor(uPlus);
- if (_patch1->_PI->periodic[1] == 1)
- vPlus -= floor(vPlus);
- //Msg::Info(" F : %g,%g,%g",F[0],F[1],F[2]);
- //Msg::Info(" R : %g,%g,%g",x0,y0,z0);
- double x0rPlus, y0rPlus, z0rPlus;
- _patch0->F(rPlus,v0,x0rPlus,y0rPlus,z0rPlus);
- double x1uPlus, y1uPlus, z1uPlus;
- _patch1->F(uPlus,v,x1uPlus,y1uPlus,z1uPlus);
- double x1vPlus, y1vPlus, z1vPlus;
- _patch1->F(u,vPlus,x1vPlus,y1vPlus,z1vPlus);
+ }
+ else
+ v1 = _SP1[1] + t * (_EP1[1] - _SP1[1]);
- double Df[3][3];
- Df[0][0] = - (x0rPlus - x0) / _h;
- Df[0][1] = (x1uPlus - x1) / _h;
- Df[0][2] = (x1vPlus - x1) / _h;
- Df[1][0] = - (y0rPlus - y0) / _h;
- Df[1][1] = (y1uPlus - y1) / _h;
- Df[1][2] = (y1vPlus - y1) / _h;
- Df[2][0] = - (z0rPlus - z0) / _h;
- Df[2][1] = (z1uPlus - z1) / _h;
- Df[2][2] = (z1vPlus - z1) / _h;
-
- //Msg::Info("D1: %g,%g,%g",Df[0][0],Df[0][1],Df[0][2]);
- //Msg::Info("D2: %g,%g,%g",Df[1][0],Df[1][1],Df[1][2]);
- //Msg::Info("D3: %g,%g,%g",Df[2][0],Df[2][1],Df[2][2]);
+ //printf("%d,%d,%g,%g: %g,%g,%g,%g\n",_patch1->_PI->periodic[0],
+ // _patch1->_PI->periodic[1],std::abs(_SP1[0]-_EP1[0]),
+ // std::abs(_SP1[1]-_EP1[1]),u0,v0,u1,v1);
- double det =
- Df[0][0] * (Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1]) +
- Df[0][1] * (Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2]) +
- Df[0][2] * (Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0]);
+ double x0, y0, z0;
+ _patch0->F(u0,v0,x0,y0,z0);
- //Msg::Info("det: %g",det);
-
- double update[3];
- update[0] =
- (Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1]) * F[0] +
- (Df[0][2] * Df[2][1] - Df[0][1] * Df[2][2]) * F[1] +
- (Df[0][1] * Df[1][2] - Df[0][2] * Df[1][1]) * F[2];
- update[1] =
- (Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2]) * F[0] +
- (Df[0][0] * Df[2][2] - Df[0][2] * Df[2][0]) * F[1] +
- (Df[0][2] * Df[1][0] - Df[0][0] * Df[1][2]) * F[2];
- update[2] =
- (Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0]) * F[0] +
- (Df[0][1] * Df[2][0] - Df[0][0] * Df[2][1]) * F[1] +
- (Df[0][0] * Df[1][1] - Df[0][1] * Df[1][0]) * F[2];
-
- //Msg::Info("U: %g,%g,%g",update[0],update[1],update[2]);
-
- r -= update[0] / det;
- u -= update[1] / det;
- v -= update[2] / det;
-
- //Msg::Info("U/det: %g,%g,%g",update[0]/det,update[1]/det,update[2]/det);
-
- //Msg::Info("uv: %g,%g,%g",r,u,v);
-
- if (_patch0->_PI->periodic[0] == 1)
- r -= std::floor(r);
- if (_patch1->_PI->periodic[0] == 1)
- u -= std::floor(u);
- if (_patch1->_PI->periodic[1] == 1)
- v -= floor(v);
-
- //Msg::Info("UV : %g,%g,%g",r,u,v);
-
- //rPlus = r + _h;
- //uPlus = u + _h;
- //vPlus = v + _h;
+ //Msg::Info("%g,%g,%g",x0,y0,z0);
- _patch0->F(r,v0,x0,y0,z0);
- _patch1->F(u,v,x1,y1,z1);
-
+ double x1, y1, z1;
+ _patch1->F(u1,v1,x1,y1,z1);
+
+ //Msg::Info("%g,%g,%g",x1,y1,z1);
+
+ double r,u,v,rPlus,uPlus,vPlus;
+
+ r = u0;
+ u = u1;
+ v = v1;
+
+ double F[3];
F[0] = x1 - x0;
F[1] = y1 - y0;
F[2] = z1 - z0;
+
+ while ((std::abs(F[0])>_tol) || (std::abs(F[1])>_tol) ||
+ (std::abs(F[2])>_tol)) {
+ rPlus = r + _h;
+ uPlus = u + _h;
+ vPlus = v + _h;
+ if (_patch0->_PI->periodic[0] == 1)
+ rPlus -= std::floor(rPlus);
+ if (_patch1->_PI->periodic[0] == 1)
+ uPlus -= std::floor(uPlus);
+ if (_patch1->_PI->periodic[1] == 1)
+ vPlus -= floor(vPlus);
+ //Msg::Info(" F : %g,%g,%g",F[0],F[1],F[2]);
+ //Msg::Info(" R : %g,%g,%g",x0,y0,z0);
+ double x0rPlus, y0rPlus, z0rPlus;
+ _patch0->F(rPlus,v0,x0rPlus,y0rPlus,z0rPlus);
+ double x1uPlus, y1uPlus, z1uPlus;
+ _patch1->F(uPlus,v,x1uPlus,y1uPlus,z1uPlus);
+ double x1vPlus, y1vPlus, z1vPlus;
+ _patch1->F(u,vPlus,x1vPlus,y1vPlus,z1vPlus);
+
+ double Df[3][3];
+ Df[0][0] = - (x0rPlus - x0) / _h;
+ Df[0][1] = (x1uPlus - x1) / _h;
+ Df[0][2] = (x1vPlus - x1) / _h;
+ Df[1][0] = - (y0rPlus - y0) / _h;
+ Df[1][1] = (y1uPlus - y1) / _h;
+ Df[1][2] = (y1vPlus - y1) / _h;
+ Df[2][0] = - (z0rPlus - z0) / _h;
+ Df[2][1] = (z1uPlus - z1) / _h;
+ Df[2][2] = (z1vPlus - z1) / _h;
+
+ //Msg::Info("D1: %g,%g,%g",Df[0][0],Df[0][1],Df[0][2]);
+ //Msg::Info("D2: %g,%g,%g",Df[1][0],Df[1][1],Df[1][2]);
+ //Msg::Info("D3: %g,%g,%g",Df[2][0],Df[2][1],Df[2][2]);
+
+ double det =
+ Df[0][0] * (Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1]) +
+ Df[0][1] * (Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2]) +
+ Df[0][2] * (Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0]);
+
+ //Msg::Info("det: %g",det);
+
+ double update[3];
+ update[0] =
+ (Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1]) * F[0] +
+ (Df[0][2] * Df[2][1] - Df[0][1] * Df[2][2]) * F[1] +
+ (Df[0][1] * Df[1][2] - Df[0][2] * Df[1][1]) * F[2];
+ update[1] =
+ (Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2]) * F[0] +
+ (Df[0][0] * Df[2][2] - Df[0][2] * Df[2][0]) * F[1] +
+ (Df[0][2] * Df[1][0] - Df[0][0] * Df[1][2]) * F[2];
+ update[2] =
+ (Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0]) * F[0] +
+ (Df[0][1] * Df[2][0] - Df[0][0] * Df[2][1]) * F[1] +
+ (Df[0][0] * Df[1][1] - Df[0][1] * Df[1][0]) * F[2];
+
+ //Msg::Info("U: %g,%g,%g",update[0],update[1],update[2]);
+
+ r -= update[0] / det;
+ u -= update[1] / det;
+ v -= update[2] / det;
+
+ //Msg::Info("U/det: %g,%g,%g",update[0]/det,update[1]/det,update[2]/det);
+
+ //Msg::Info("uv: %g,%g,%g",r,u,v);
+
+ if (_patch0->_PI->periodic[0] == 1)
+ r -= std::floor(r);
+ if (_patch1->_PI->periodic[0] == 1)
+ u -= std::floor(u);
+ if (_patch1->_PI->periodic[1] == 1)
+ v -= floor(v);
+
+ //Msg::Info("UV : %g,%g,%g",r,u,v);
+
+ //rPlus = r + _h;
+ //uPlus = u + _h;
+ //vPlus = v + _h;
+
+ _patch0->F(r,v0,x0,y0,z0);
+ _patch1->F(u,v,x1,y1,z1);
+
+ F[0] = x1 - x0;
+ F[1] = y1 - y0;
+ F[2] = z1 - z0;
+ }
+ x = x0; y = y0; z = z0;
+ }
+ else if (_patch0) {
+ if (_II->edgeNumber == 0)
+ _patch0->F(t,0.,x,y,z);
+ else if (_II->edgeNumber == 1)
+ _patch0->F(1.,t,x,y,z);
+ else if (_II->edgeNumber == 2)
+ _patch0->F(1.-t,1.,x,y,z);
+ else if (_II->edgeNumber == 3)
+ _patch0->F(0.,1.-t,x,y,z);
+ else
+ Msg::Error("Unknown edge number.");
}
- x = x0; y = y0; z = z0;
+ else if (_patch1) {
+ if (_II->edgeNumber == 0)
+ _patch1->F(t,0.,x,y,z);
+ else if (_II->edgeNumber == 1)
+ _patch1->F(1.,t,x,y,z);
+ else if (_II->edgeNumber == 2)
+ _patch1->F(1.-t,1.,x,y,z);
+ else if (_II->edgeNumber == 3)
+ _patch1->F(0.,1.-t,x,y,z);
+ }
+ else
+ Msg::Error("No patches to find the intersection curve.");
}
bool Curve::Inverse(double x,double y,double z,double &t)
{
- double u0, v0;
- _patch0->Inverse(x,y,z,u0,v0);
-
- if ((std::abs(_SP0[1]-_EP0[1])<1.e-12) &&
- (_patch0->_PI->periodic[1] == 1)) {
- t = (v0 - _SP0[1] > 0 ? v0 - _SP0[1] : 1. + v0 - _SP0[1]);
+ if ((_patch0) && (_patch1)) {
+ double u0, v0;
+ _patch0->Inverse(x,y,z,u0,v0);
+
+ if ((std::abs(_SP0[1]-_EP0[1])<1.e-12) &&
+ (_patch0->_PI->periodic[1] == 1)) {
+ t = (v0 - _SP0[1] > 0 ? v0 - _SP0[1] : 1. + v0 - _SP0[1]);
+ }
+ else
+ t = (v0 - _SP0[1]) / (_EP0[1] - _SP0[1]);
+ }
+ else if (_patch0) {
+ double u,v;
+ _patch0->Inverse(x,y,z,u,v);
+ if (_II->edgeNumber == 0)
+ t = u;
+ else if (_II->edgeNumber == 1)
+ t = v;
+ else if (_II->edgeNumber == 2)
+ t = 1. - u;
+ else if (_II->edgeNumber == 3)
+ t = 1. - v;
+ }
+ else if (_patch1) {
+ double u,v;
+ _patch1->Inverse(x,y,z,u,v);
+ if (_II->edgeNumber == 0)
+ t = u;
+ else if (_II->edgeNumber == 1)
+ t = v;
+ else if (_II->edgeNumber == 2)
+ t = 1. - u;
+ else if (_II->edgeNumber == 3)
+ t = 1. - v;
}
else
- t = (v0 - _SP0[1]) / (_EP0[1] - _SP0[1]);
+ Msg::Error("No patches to find the intersection curve.");
if ((t < 0.) || (t > 1.))
return false;
diff --git a/contrib/FourierModel/Curve.h b/contrib/FourierModel/Curve.h
index 8606061317..f1934bcd53 100644
--- a/contrib/FourierModel/Curve.h
+++ b/contrib/FourierModel/Curve.h
@@ -7,6 +7,7 @@
// The base class for the patches
class Curve {
private:
+ int _edgeNum;
double _h, _tol;
protected:
// Patches
diff --git a/contrib/FourierModel/FM_Edge.cpp b/contrib/FourierModel/FM_Edge.cpp
index eb87fa4d53..e8ffbab5cb 100644
--- a/contrib/FourierModel/FM_Edge.cpp
+++ b/contrib/FourierModel/FM_Edge.cpp
@@ -76,19 +76,19 @@ void FM_Edge::Dfdt(double t, double &x, double &y, double &z)
}
else
tRescaled = tStart + t * (tEnd - tStart);
- if (t+0.5*h > 1.) {
+ if (tRescaled+0.5*h > tEnd) {
_curve->F(tRescaled, xPlus, yPlus, zPlus);
- double tMinus = tStart + (t - h) * (tEnd - tStart);
+ double tMinus = tRescaled - h;
_curve->F(tMinus, xMinus, yMinus, zMinus);
}
- else if (t-0.5*h < 0.) {
+ else if (tRescaled-0.5*h < tStart) {
_curve->F(tRescaled, xMinus, yMinus, zMinus);
- double tPlus = tStart + (t + h) * (tEnd - tStart);
+ double tPlus = tRescaled + h;
_curve->F(tPlus, xPlus, yPlus, zPlus);
}
else {
- double tPlus = tStart + (t + 0.5 * h) * (tEnd - tStart);
- double tMinus = tStart + (t - 0.5 * h) * (tEnd - tStart);
+ double tPlus = tRescaled + 0.5 * h;
+ double tMinus = tRescaled - 0.5 * h;
_curve->F(tPlus, xPlus, yPlus, zPlus);
_curve->F(tMinus, xMinus, yMinus, zMinus);
}
@@ -121,19 +121,19 @@ void FM_Edge::Dfdfdtdt(double t, double &x, double &y, double &z)
}
else
tRescaled = tStart + t * (tEnd - tStart);
- if (t+0.5*h > 1.) {
+ if (tRescaled+0.5*h > tEnd) {
Dfdt(tRescaled, xPlus, yPlus, zPlus);
- double tMinus = tStart + (t - h) * (tEnd - tStart);
+ double tMinus = tRescaled - h;
Dfdt(tMinus, xMinus, yMinus, zMinus);
}
- else if (t-0.5*h < 0.) {
+ else if (tRescaled-0.5*h < tStart) {
Dfdt(tRescaled, xMinus, yMinus, zMinus);
- double tPlus = tStart + (t + h) * (tEnd - tStart);
+ double tPlus = tRescaled + h;
Dfdt(tPlus, xPlus, yPlus, zPlus);
}
else {
- double tPlus = tStart + (t + 0.5 * h) * (tEnd - tStart);
- double tMinus = tStart + (t - 0.5 * h) * (tEnd - tStart);
+ double tPlus = tRescaled + 0.5 * h;
+ double tMinus = tRescaled - 0.5 * h;
Dfdt(tPlus, xPlus, yPlus, zPlus);
Dfdt(tMinus, xMinus, yMinus, zMinus);
}
diff --git a/contrib/FourierModel/FM_Info.h b/contrib/FourierModel/FM_Info.h
index 4eb3a0071f..9c723b6a2b 100644
--- a/contrib/FourierModel/FM_Info.h
+++ b/contrib/FourierModel/FM_Info.h
@@ -29,8 +29,9 @@ class IntersectionInfo {
struct {
int patchTag;
} intersectingPatches[2];
+ int edgeNumber;
- IntersectionInfo() {}
+ IntersectionInfo() : tag(-1), edgeNumber(-1) {}
virtual ~IntersectionInfo() {}
};
diff --git a/contrib/FourierModel/FM_Reader.cpp b/contrib/FourierModel/FM_Reader.cpp
index 7f5398857d..994bca5137 100644
--- a/contrib/FourierModel/FM_Reader.cpp
+++ b/contrib/FourierModel/FM_Reader.cpp
@@ -42,6 +42,7 @@ FM_Reader::FM_Reader(const char* fn)
}
}
}
+
for (int i=0;i<_nIntersections;i++) {
_intersectionList.push_back(new IntersectionInfo);
InputFile >> _intersectionList[i]->tag;
@@ -51,10 +52,13 @@ FM_Reader::FM_Reader(const char* fn)
>> _intersectionList[i]->EP[2];
InputFile >> _intersectionList[i]->intersectingPatches[0].patchTag;
InputFile >> _intersectionList[i]->intersectingPatches[1].patchTag;
+ if ((_intersectionList[i]->intersectingPatches[0].patchTag < 0) ||
+ (_intersectionList[i]->intersectingPatches[1].patchTag < 0))
+ InputFile >> _intersectionList[i]->edgeNumber;
}
_patch.resize(_nPatches);
- _intersection.resize(_nIntersections);
+ _intersection.resize(_nIntersections,0);
for (int i=0;i<_nPatches;i++) {
PatchInfo* PI = _patchList[i];
@@ -109,7 +113,12 @@ Patch* FM_Reader::GetPatch(int tag)
Curve* FM_Reader::GetIntersection(int tag)
{
- for (int i=0;i<_intersection.size();i++)
- if (_intersection[i]->GetTag() == tag)
- return _intersection[i];
+ Curve* curve = 0;
+ for (int i=0;i<_intersection.size();i++) {
+ if (_intersection[i]->GetTag() == tag) {
+ curve = _intersection[i];
+ break;
+ }
+ }
+ return curve;
}
--
GitLab