diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 6a639cde7068a731d86ea33d4bf93a42cd0ca420..591228d545322973170a92c1e5ee1b62e39210ad 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -1,3 +1,4 @@
+
// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
@@ -26,6 +27,118 @@
#include "convexCombinationTerm.h" // #
#include "BasisFactory.h"
+static inline void functionShapes(int p, double Xi[2], double* phi)
+{
+
+ switch(p){
+
+ case 1:
+ phi[0] = 1. - Xi[0] - Xi[1];
+ phi[1] = Xi[0];
+ phi[2] = Xi[1];
+ break;
+
+ case 2:
+ phi[0] = 1. - 3.*(Xi[0]+Xi[1]) + 2.*(Xi[0]+Xi[1])*(Xi[0]+Xi[1]);
+ phi[1] = Xi[0]*(2.*Xi[0]-1);
+ phi[2] = Xi[1]*(2.*Xi[1]-1);
+ phi[3] = 4.*Xi[0]*(1.-Xi[0]-Xi[1]);
+ phi[4] = 4.*Xi[0]*Xi[1];
+ phi[5] = 4.*Xi[1]*(1.-Xi[0]-Xi[1]);
+ break;
+
+ default:
+ Msg::Error("(discreteDiskFace) static inline functionShapes, only first and second order available; order %d requested.",p);
+ break;
+ }
+
+}
+
+
+static inline void derivativeShapes(int p, double Xi[2], double phi[6][2])
+{
+
+ switch(p){
+
+ case 1:
+ phi[0][0] = -1. ; phi[0][1] = -1.;
+ phi[1][0] = 1. ; phi[1][1] = 0.;
+ phi[2][0] = 0. ; phi[2][1] = 1.;
+
+ break;
+
+ case 2:
+ phi[0][0] = -3. + 4.*(Xi[0]+Xi[1]) ; phi[0][1] = -3. + 4.*(Xi[0]+Xi[1]);
+ phi[1][0] = 4.*Xi[0] - 1. ; phi[1][1] = 0.;
+ phi[2][0] = 0. ; phi[2][1] = 4.*Xi[1] - 1.;
+ phi[3][0] = 4. - 8.*Xi[0] - 4.*Xi[1] ; phi[3][1] = -4.*Xi[0];
+ phi[4][0] = 4.*Xi[1] ; phi[4][1] = 4.*Xi[0];
+ phi[5][0] = -4.*Xi[1] ; phi[5][1] = 4. - 4.*Xi[0] - 8.*Xi[1];
+ break;
+
+ default:
+ Msg::Error("(discreteDiskFace) static inline derivativeShapes, only first and second order available; order %d requested.",p);
+ break;
+ }
+
+}
+
+static inline void uv2xi(discreteDiskFaceTriangle* my_ddft, double U[2], double Xi[2]){
+
+ double M[2][2], R[2];
+ const SPoint3 p0 = my_ddft->p[0];
+ const SPoint3 p1 = my_ddft->p[1];
+ const SPoint3 p2 = my_ddft->p[2];
+ M[0][0] = p1.x() - p0.x();
+ M[0][1] = p2.x() - p0.x();
+ M[1][0] = p1.y() - p0.y();
+ M[1][1] = p2.y() - p0.y();
+ R[0] = (U[0] - p0.x());
+ R[1] = (U[1] - p0.y());
+ sys2x2(M, R, Xi);
+
+ if (my_ddft->tri->getPolynomialOrder() == 2){
+
+ int iter = 1, maxiter = 20;
+ double error = 1., tol = 1.e-6;
+ while (error > tol && iter < maxiter){// Newton-Raphson
+
+ double fs[6];// phi_i
+ functionShapes(2,Xi,fs);
+ double ds[6][2];// dPhi_i/dXi_j
+ derivativeShapes(2,Xi,ds);
+
+ double un[2] = {0.,0.};
+ double jac[2][2] = {{0.,0.},{0.,0.}}; // d(u,v)/d(xi,eta)
+ for (int i=0; i<6; i++){
+ double ui[2] = {my_ddft->p[i].x(),my_ddft->p[i].y()};
+ for(int k=0; k<2; k++){
+ un[k] += ui[k] * fs[i];
+ for (int j=0; j<2; j++)
+ jac[k][j] += ui[k] * ds[i][j];
+ }
+ }
+
+ double inv[2][2];
+ inv2x2(jac,inv);
+ double xin[2] = {Xi[0],Xi[1]};
+ error = 0.;
+ for (int i=0; i<2; i++){
+ for (int j=0; j<2; j++)
+ xin[i] += inv[i][j] * (U[j] - un[j]);
+ error += (xin[i] - Xi[i])*(xin[i] - Xi[i]);
+ }
+
+ error = sqrt(error);
+ Xi[0] = xin[0];
+ Xi[1] = xin[1];
+ iter++;
+
+ } // end Newton-Raphson
+ }// end order 2
+}
+
+
// The three things that are mandatory to manipulate octrees (octree in (u;v)).
static void discreteDiskFaceBB(void *a, double*mmin, double*mmax)
{ // called once by buildOct()
@@ -63,9 +176,11 @@ static int discreteDiskFaceInEle(void *a, double*c)// # mark
{ // called once by buildOct()
discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a;
-
- double X[2];
- const double eps = 1.e-8;
+ double Xi[2];
+ double U[2] = {c[0],c[1]};
+ uv2xi(t,U,Xi);
+ double eps = 1e-8;
+ /*
if (t->tri->getPolynomialOrder() == 1){
@@ -95,8 +210,8 @@ static int discreteDiskFaceInEle(void *a, double*c)// # mark
t6.xyz2uvw(xyz,X);
}
-
- if(X[0] > -eps && X[1] > -eps && 1. - X[0] - X[1] > -eps)
+ */
+ if(Xi[0] > -eps && Xi[1] > -eps && 1. - Xi[0] - Xi[1] > -eps)
return 1;
return 0;
@@ -146,7 +261,7 @@ discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh, int
tagElementOrder = -1;
Msg::Error("discreteDiskFace:: only p=1 or p=2 allowed");
}
- mynodalbasis = BasisFactory::getNodalBasis(tagElementOrder);
+ //mynodalbasis = BasisFactory::getNodalBasis(tagElementOrder);
std::map<MVertex*,MVertex*> v2v;// mesh vertex |-> face vertex
std::map<MEdge,MVertex*,Less_Edge> ed2nodes; // edge to interior node(s)
@@ -477,11 +592,18 @@ bool discreteDiskFace::parametrize() const
void discreteDiskFace::getTriangleUV(const double u,const double v,
discreteDiskFaceTriangle **mt,
- double &_u, double &_v)const{
+ double &_xi, double &_eta)const{
double uv[3] = {u,v,0.};
*mt = (discreteDiskFaceTriangle*) Octree_Search(uv,oct);
if (!(*mt)) return;
+
+ double Xi[2];
+ double U[2] = {u,v};
+ uv2xi(*mt,U,Xi);
+ _xi = Xi[0];
+ _eta = Xi[1];
+ /*
if (_order == 1){
double M[2][2],X[2],R[2];
const SPoint3 p0 = (*mt)->p[0];
@@ -518,7 +640,8 @@ void discreteDiskFace::getTriangleUV(const double u,const double v,
_u = uv[0];// xi
_v = uv[1];// eta
- }
+ }
+ */
}
void discreteDiskFace::checkOrientationUV(){
@@ -549,19 +672,21 @@ void discreteDiskFace::checkOrientationUV(){
// (u;v) |-> < (x;y;z); GFace; (u;v) >
GPoint discreteDiskFace::point(const double par1, const double par2) const
{
- double U,V;
+ double xi,eta;
double par[2] = {par1,par2};
discreteDiskFaceTriangle* dt;
- getTriangleUV(par1,par2,&dt,U,V);// return Xi,Eta
+ getTriangleUV(par1,par2,&dt,xi,eta);// return Xi,Eta
+ double Xi[2] = {xi,eta};
if (!dt) {
GPoint gp = GPoint(1.e22,1.e22,1.e22,_parent,par);
gp.setNoSuccess();
return gp;
}
-
+
//polynomialBasis myPolynomial = polynomialBasis(dt->tri->getTypeForMSH());
double eval[_N];
- mynodalbasis->f(U,V,0.,eval);//#mark
+ //mynodalbasis->f(U,V,0.,eval);//#mark
+ functionShapes(_order,Xi,eval);
double X=0,Y=0,Z=0;
for(int io=0; io<_N; io++){
X += dt->tri->getVertex(io)->x()*eval[io];
@@ -626,9 +751,9 @@ double discreteDiskFace::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVec
Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 ¶m) const
{
- double U,V;
+ double xi,eta;
discreteDiskFaceTriangle* ddft;
- getTriangleUV(param.x(),param.y(),&ddft,U,V);
+ getTriangleUV(param.x(),param.y(),&ddft,xi,eta);
MTriangle* tri = NULL;
if (ddft) tri = ddft->tri;
@@ -637,10 +762,10 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 ¶m) const
return Pair<SVector3, SVector3>(SVector3(1,0,0), SVector3(0,1,0));
}
-
- double df[_N][3];
- mynodalbasis->df(U,V,0.,df);
-
+ double Xi[2] = {xi,eta};
+ double df[_N][2];
+ //mynodalbasis->df(U,V,0.,df);
+ derivativeShapes(_order,Xi,df);
double dxdxi[3][2] = {{0.,0.},{0.,0.},{0.,0.}};
double dudxi[2][2] = {{0.,0.},{0.,0.}};
@@ -885,6 +1010,4 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
}
-
-
#endif
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index 499dcf875c3b77360ab56a7d08bde15db771e8a0..61a845d61d923a3ab636274e20ada61a9588b8dc 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -96,7 +96,7 @@ class discreteDiskFace : public GFace {
mutable discreteDiskFaceTriangle *_ddft;
mutable Octree *oct;
mutable std::vector<double> _coords;
- const nodalBasis* mynodalbasis;
+ // const nodalBasis* mynodalbasis;
};
#endif