diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index c5f62c5bfc68cda651342c895fcef362c5733be8..47085b4491741e7547531a8e8fa91f587eb5d2d6 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -3,6 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+#include <queue> // #mark
#include <stdlib.h>
#include "GmshConfig.h"
@@ -23,18 +24,33 @@
#include "laplaceTerm.h"
#include "convexLaplaceTerm.h"
#include "convexCombinationTerm.h" // #
+#include "BasisFactory.h"
// 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()
discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a;
- mmin[0] = std::min(std::min(t->p1.x(), t->p2.x()), t->p3.x());
- mmin[1] = std::min(std::min(t->p1.y(), t->p2.y()), t->p3.y());
- mmax[0] = std::max(std::max(t->p1.x(), t->p2.x()), t->p3.x());
- mmax[1] = std::max(std::max(t->p1.y(), t->p2.y()), t->p3.y());
+
+
+ mmin[0] = std::min(std::min(t->p[0].x(), t->p[1].x()), t->p[2].x());
+ mmin[1] = std::min(std::min(t->p[0].y(), t->p[1].y()), t->p[2].y());
+ mmax[0] = std::max(std::max(t->p[0].x(), t->p[1].x()), t->p[2].x());
+ mmax[1] = std::max(std::max(t->p[0].y(), t->p[1].y()), t->p[2].y());
+
+
+ if (t->tri->getPolynomialOrder() == 2){
+ for (int i=3; i<6; i++){
+ int j = i==5 ? 0 : i-2;
+ double du = t->p[i].x() - (t->p[i-3].x() + t->p[j].x())/2.;
+ double dv = t->p[i].y() - (t->p[i-3].y() + t->p[j].y())/2.;
+ mmin[0] = std::min(mmin[0],t->p[i].x()+du);
+ mmin[1] = std::min(mmin[1],t->p[i].y()+dv);
+ mmax[0] = std::max(mmax[0],t->p[i].x()+du);
+ mmax[1] = std::max(mmax[1],t->p[i].y()+dv);
+ }
+ }
mmin[2] = -1;
mmax[2] = 1;
-
const double dx = mmax[0] - mmin[0];
const double dy = mmax[1] - mmin[1];
const double eps = 0.0;//1.e-12;
@@ -47,26 +63,44 @@ static void discreteDiskFaceBB(void *a, double*mmin, double*mmax)
static void discreteDiskFaceCentroid(void *a, double*c)
{ // called once by buildOct()
discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a;
- c[0] = (t->p1.x() + t->p2.x() + t->p3.x()) / 3.0;
- c[1] = (t->p1.y() + t->p2.y() + t->p3.y()) / 3.0;
+ c[0] = (t->p[0].x() + t->p[1].x() + t->p[2].x()) / 3.0;
+ c[1] = (t->p[0].y() + t->p[1].y() + t->p[2].y()) / 3.0;
c[2] = 0.0;
}
-static int discreteDiskFaceInEle(void *a, double*c)
+static int discreteDiskFaceInEle(void *a, double*c)// # mark
{ // called once by buildOct()
+
discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a;
- double M[2][2], R[2], X[2];
+
+ double X[2];
const double eps = 1.e-8;
- const SPoint3 p0 = t->p1;
- const SPoint3 p1 = t->p2;
- const SPoint3 p2 = t->p3;
- 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] = (c[0] - p0.x());
- R[1] = (c[1] - p0.y());
- sys2x2(M, R, X);
+
+
+
+ if (t->tri->getPolynomialOrder() == 1){
+
+ double M[2][2], R[2];
+ const SPoint3 p0 = t->p[0];
+ const SPoint3 p1 = t->p[1];
+ const SPoint3 p2 = t->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] = (c[0] - p0.x());
+ R[1] = (c[1] - p0.y());
+ sys2x2(M, R, X);
+ }
+ else{
+
+ double uv[3] = {c[0],c[1],0.};
+ double Xi[3];
+ t->tri->xyz2uvw(uv,Xi);
+ X[0] = Xi[0];
+ X[1] = Xi[1];
+ }
+
if(X[0] > -eps && X[1] > -eps && 1. - X[0] - X[1] > -eps)
return 1;
@@ -99,97 +133,301 @@ static bool orderVertices(const double &tot_length, const std::vector<MVertex*>
}
/*BUILDER*/
-discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh) :
+discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh, int p) :
GFace(gf->model(),123), _parent (gf)
{
- if (mesh.empty())return;
- std::map<MVertex*,MVertex*> v2v;
- for (unsigned int i=0;i<mesh.size();i++){
- MVertex *vs[3] = {NULL, NULL, NULL};
- for (int j=0;j<3;j++){
- MVertex *v = mesh[i]->getVertex(j);
+
+ _order = p;
+ _N = (p+1)*(p+2)/2;
+
+ int tagElementOrder;
+ switch (_order){
+ case 1:
+ tagElementOrder = MSH_TRI_3;
+ break;
+ case 2:
+ tagElementOrder = MSH_TRI_6;
+ break;
+ case 3:
+ tagElementOrder = MSH_TRI_10;
+ break;
+ default:
+ Msg::Error("discreteDiskFace:: builder, tag is not (yet) covered for this order :-) ");
+ }
+
+ mynodalbasis = BasisFactory::getNodalBasis(tagElementOrder);
+ fullMatrix<double> myNodes;
+ mynodalbasis->getReferenceNodes(myNodes);
+
+ mybezierbasis = BasisFactory::getBezierBasis(TYPE_TRI,_order);
+
+ std::map<MVertex*,MVertex*> v2v;// mesh vertex |-> face vertex
+ std::map<MEdge,std::vector<MVertex*>,Less_Edge> ed2nodes; // edge to interior node(s)
+ for (unsigned int i=0;i<mesh.size();i++){ // triangle by triangle
+
+ std::vector<MVertex*> vs; // MTriangle vertices
+ for (unsigned int j=0; j<3; j++){ // loop over vertices AND edges of the current triangle
+
+ MVertex *v = mesh[i]->getVertex(j);// firstly, edge vertices
if (v->onWhat() == gf) {
std::map<MVertex*,MVertex*> :: iterator it = v2v.find(v);
if (it == v2v.end()){
MFaceVertex *vv = new MFaceVertex ( v->x(), v->y(), v->z(), this, 0, 0);
v2v[v] = vv;
discrete_vertices.push_back(vv);
- vs[j] = vv;
+ vs.push_back(vv);
}
- else vs[j] = it->second;
+ else vs.push_back(it->second);
}
else if (v->onWhat()->dim() == 1){
if (v->onWhat()->geomType() == DiscreteCurve){
discreteEdge *de = dynamic_cast<discreteEdge*> (v->onWhat());
- vs[j] = de->getGeometricalVertex(v);
+ vs.push_back(de->getGeometricalVertex(v));
}
+
else Msg::Fatal("fatality");
}
- else vs[j] = v;
- }
- discrete_triangles.push_back(new MTriangle (vs[0], vs[1], vs[2]));
- }
+ else vs.push_back(v);
+
+
+ }// end loop over vertices and edges
+ if (_order >= 2) {// then, interior nodes :-)
+ double Xi, Eta, X, Y,Z;
+
+ for (unsigned int ie=0; ie<3; ie++){ // firstly, edge interior nodes :-)
+
+ MEdge me = mesh[i]->getEdge(ie); // check if edge already visited
+ if(ed2nodes.find(me) == ed2nodes.end()){
+
+ for(unsigned int iev=3+ie*(p-1);iev<3+(ie+1)*(p-1);iev++){
+ Xi = myNodes(iev,0);
+ Eta = myNodes(iev,1);
+
+ X = vs[0]->x()*(1.-Xi-Eta) + vs[1]->x()*Xi + vs[2]->x()*Eta;
+ Y = vs[0]->y()*(1.-Xi-Eta) + vs[1]->y()*Xi + vs[2]->y()*Eta;
+ Z = vs[0]->z()*(1.-Xi-Eta) + vs[1]->z()*Xi + vs[2]->z()*Eta;
+
+ MVertex* mv = new MVertex(X,Y,Z,gf);
+ vs.push_back(mv);
+ ed2nodes[me].push_back(mv);
+ }
+ }
+ else{
+ std::vector<MVertex*> vecmv = ed2nodes[me];
+ for(unsigned int iev=0; iev<(p-1); iev++)
+ vs.push_back(vecmv[iev]);
+ }
+ }
+
+ for (unsigned int iv=3*_order; iv<_N; iv++){ // then, volume nodes
+
+ Xi = myNodes(iv,0);
+ Eta = myNodes(iv,1);
+
+ X = vs[0]->x()*(1.-Xi-Eta) + vs[1]->x()*Xi + vs[2]->x()*Eta;
+ Y = vs[0]->y()*(1.-Xi-Eta) + vs[1]->y()*Xi + vs[2]->y()*Eta;
+ Z = vs[0]->z()*(1.-Xi-Eta) + vs[1]->z()*Xi + vs[2]->z()*Eta;
+
+ vs.push_back(new MVertex(X,Y,Z,gf));
+ }
+ }// end order > 2
+
+ if (_order==1)
+ discrete_triangles.push_back(new MTriangle (vs));
+
+ if (_order==2)
+ discrete_triangles.push_back(new MTriangle6 (vs));
+
+ if (_order > 2)
+ discrete_triangles.push_back(new MTriangleN (vs,(char)_order));
+
+ for(unsigned int check=0; check<_N; check++)
+ std::cout << "(" << vs[check]->x() << ";" << vs[check]->y() << ")\t";
+ std::cout << "\n";
+ for(unsigned int check=0; check<_N; check++)
+ std::cout << "(" << myNodes(check,0) << ";" << myNodes(check,1) << ")\t";
+ std::cout << "\n";
+
+ }// end loop over triangles
+
+ // delete MTriangle* and MVertex ? Y E S # mark
+
- // triangles = mesh;
uv_kdtree = NULL;
- kdtree = NULL;
- // checkOrientation(); // mark, need to check after parametrization ? --> no ! (3D, cf. Seb)
+ kdtree = NULL;
+ std::cout << "buildAllNodes" << std::endl;
buildAllNodes();
+ std::cout << "getBoundingEdges" << std::endl;
getBoundingEdges();
+ std::cout << "orderVertices" << std::endl;
orderVertices(_totLength, _U0, _coords);
+ std::cout << "parametrize" << std::endl;
parametrize();
+ std::cout << "buildOct" << std::endl;
buildOct();
- // putOnView();
+ std::cout << "checkOrientation" << std::endl;
+ checkOrientationUV();
+ std::cout << "putOnView" << std::endl;
+ putOnView();
+ std::cout << "discreteDiskFace: The parametrization of triangulation " << gf->tag() << " is ending." << std::endl;
}
+
void discreteDiskFace::getBoundingEdges()
{
+
+ // first of all, all the triangles have to be oriented in the same way
+ std::map<MEdge,std::vector<MElement*>,Less_Edge> ed2tri; // edge to 1 or 2 triangle(s)
+
+ for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
+ MElement *e = discrete_triangles[i];
+ for(int j = 0; j < e->getNumEdges() ; j++){
+ MEdge ed = e->getEdge(j);
+ ed2tri[ed].push_back(e);
+ }
+ }
+
+ // element to its neighbors
+ std::map<MElement*,std::vector<MElement*> > neighbors;
+ for (unsigned int i=0; i<discrete_triangles.size(); ++i){
+ MElement* e = discrete_triangles[i];
+ for( unsigned int j=0; j<e->getNumEdges(); j++){ // #fixme: efficiency could be improved by setting neighbors mutually
+ std::vector<MElement*> my_mt = ed2tri[e->getEdge(j)];
+ if (my_mt.size() > 1){// my_mt.size() = {1;2}
+ MElement* neighTri = my_mt[0] == e ? my_mt[1] : my_mt[0];
+ neighbors[e].push_back(neighTri);
+ }
+ }
+ }
+
+
+ // queue: first in, first out
+ std::queue<MElement*> checkList; // element for reference orientation
+ std::queue< std::vector<MElement*> > checkLists; // corresponding neighbor element to be checked for its orientation
+ std::queue<MElement*> my_todo; // todo list
+ std::map<MElement*,bool> check_todo; // help to complete todo list
+
+ my_todo.push(discrete_triangles[0]);
+ check_todo[discrete_triangles[0]]=true;
+ while(!my_todo.empty()){
+
+ MElement* myMT = my_todo.front();
+ my_todo.pop();
+
+ std::vector<MElement*> myV = neighbors[myMT];
+ std::vector<MElement*> myInsertion;
+
+ checkList.push(myMT);
+
+ for(unsigned int i=0; i<myV.size(); ++i){
+ if (check_todo.find(myV[i]) == check_todo.end()){
+ myInsertion.push_back(myV[i]);
+ check_todo[myV[i]] = true;
+ my_todo.push(myV[i]);
+ }
+ }
+ checkLists.push(myInsertion);
+ }// end while
+
+
+ while(!checkList.empty() && !checkLists.empty()){
+ MElement* current = checkList.front();
+ checkList.pop();
+ std::vector<MElement*> neigs = checkLists.front();
+ checkLists.pop();
+ for (unsigned int i=0; i<neigs.size(); i++){
+ bool myCond = false;
+ for (unsigned int k=0; k<3; k++){
+ for (unsigned int j=0; j<3; j++){
+ if (current->getVertex(k) == neigs[i]->getVertex(j)){
+ myCond = true;
+ if (!(current->getVertex(k!=2 ?k+1 : 0 ) == neigs[i]->getVertex(j!=0 ? j-1 : 2) ||
+ current->getVertex(k!=0 ?k-1 : 2 ) == neigs[i]->getVertex(j!=2 ? j+1 : 0))){
+ neigs[i]->reverse();
+ std::cout << "discreteDiskFace: triangle " << neigs[i]->getNum() << " has been reoriented." << std::endl;
+ }
+ break;
+ }
+ }
+ if (myCond)
+ break;
+ }
+ } // end for unsigned int i
+ } // end while
+
+
+ // now, it is possible to identify the border(s) of the triangulation
// allEdges contains all edges of the boundary
std::set<MEdge,Less_Edge> allEdges;
-
+
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
MElement *e = discrete_triangles[i];
for(int j = 0; j < e->getNumEdges() ; j++){
MEdge ed = e->getEdge(j);
std::set<MEdge,Less_Edge>::iterator it = allEdges.find(ed);
- if (it == allEdges.end()) allEdges.insert(ed);
- else allEdges.erase(it);
+ if (it == allEdges.end()) allEdges.insert(ed);
+ else allEdges.erase(it);
}
}
- std::map<MVertex*,MVertex*> edges; // useless ?
- std::map<MVertex*,MVertex*> firstNode2Edge;
+ std::map<MVertex*,std::vector<MVertex*> > firstNode2Edge;
for (std::set<MEdge>::iterator ie = allEdges.begin(); ie != allEdges.end() ; ++ie) {
- MEdge ed = *ie;
- MVertex *first = ed.getVertex(0);
- MVertex *last = ed.getVertex(1);
+ MEdge ed = *ie;
- std::map<MVertex*,MVertex*>::iterator im = firstNode2Edge.find(first);
- if (im != firstNode2Edge.end()) Msg::Fatal("Incorrect topology in discreteDiskFace %d", tag());
+ std::vector<MElement *> vectri = ed2tri[ed];
+ MElement* tri = vectri[0]; // boundary edge: only one triangle :-)
+
+ std::vector<MVertex*> vecver;
+ tri->getEdgeVertices(edgeLocalNum(tri,ed),vecver);
+ MVertex *first = vecver[0];
+ MVertex *last = vecver[1];
+ vecver.erase(vecver.begin());
+ vecver.erase(vecver.begin());
- firstNode2Edge[first] = last;
+ std::map<MVertex*,std::vector<MVertex*> >::iterator im = firstNode2Edge.find(first);
+ if (im != firstNode2Edge.end()) Msg::Fatal("Incorrect topology in discreteDiskFace %d", tag());
+
+ firstNode2Edge[first] = vecver;
+ firstNode2Edge[first].push_back(last);
}
- while (!firstNode2Edge.empty()) { // quid not empty but in within the map ?
+ while (!firstNode2Edge.empty()) { // quid 'firstNode2Edge' is not empty but 'in' within the map ?
std::vector<MVertex*> loop;
double length = 0.;
- std::map<MVertex*,MVertex*>::iterator in = firstNode2Edge.begin();
+
+ std::map<MVertex*,std::vector<MVertex*> >::iterator in = firstNode2Edge.begin();
while(in != firstNode2Edge.end()) {
MVertex *first = in->first;
- MVertex *last = in->second;
- length += sqrt( (last->x()-first->x()) * (last->x()-first->x()) +
- (last->y()-first->y()) * (last->y()-first->y()) +
- (last->z()-first->z()) * (last->z()-first->z()) );
+ std::vector<MVertex*> myV = in->second;
+
loop.push_back(first);
- firstNode2Edge.erase(in);
- in = firstNode2Edge.find(last);
- }
- _loops.insert(std::make_pair(length,loop)); // it shouldn't be possible to have twice the same length ? mark
- }
+ MVertex* previous = first;
+
+ for(int i=0; i<=myV.size()-1; i++){
- // dirichlet BC
+ MVertex* current = myV[i];
+
+ std::cout << "(" << current->x() << ";" << current->y() << ")" << std::endl;
+
+ loop.push_back(current);
+
+ length += sqrt( (current->x()-previous->x()) * (current->x()-previous->x()) +
+ (current->y()-previous->y()) * (current->y()-previous->y()) +
+ (current->z()-previous->z()) * (current->z()-previous->z()) );
+
+ _loops.insert(std::make_pair(length,loop)); // it shouldn't be possible to have twice the same length ? actually, it is possible, but quite seldom #fixme ----> multimap ?
+
+ previous = current;
+
+ }
+ firstNode2Edge.erase(in);
+ in = firstNode2Edge.find(previous);
+ }// end while in
+ } // end while firstNode2Edge
+
+ // dirichlet BC
_U0 = _loops.rbegin()->second;
_totLength = _loops.rbegin()->first;
}
@@ -198,21 +436,26 @@ void discreteDiskFace::buildOct() const
{
SBoundingBox3d bb;
int count = 0;
-
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
MTriangle *t = discrete_triangles[i];
//create bounding box
- for(int j = 0; j < 3; j++){
+ for(int j = 0; j < t->getNumVertices(); j++){
std::map<MVertex*,SPoint3>::const_iterator itj = coordinates.find(t->getVertex(j));
- _coordPoints.insert(std::make_pair(t->getVertex(j)->point(), itj->second));// <tr(x;y;z);tr(u;v)>
- bb += SPoint3(itj->second.x(),itj->second.y(),0.0);
+ _coordPoints.insert(std::make_pair(t->getVertex(j)->point(), itj->second));// <tr(x;y;z);tr(u;v)>
+ if (_order == 2 && j>2){
+ int ij = j == 5 ? 0 : j-2;
+ double du = itj->second.x() - (t->getVertex(j-3)->x()+t->getVertex(ij)->x())/2.;
+ double dv = itj->second.y() - (t->getVertex(j-3)->y()+t->getVertex(ij)->y())/2.;
+ bb += SPoint3(itj->second.x()+du,itj->second.y()+dv,0.);
+ }
+ else
+ bb += SPoint3(itj->second.x(), itj->second.y(),0.);
}
count++;
}
- //ANN octree
- std::set<MVertex*> allVS; // ? useless ?
- ANNpointArray nodes = annAllocPts(count, 3);
+ // ANN octree
+ ANNpointArray nodes = annAllocPts(count, _N);
// make bounding box
SPoint3 bbmin = bb.min(), bbmax = bb.max();
@@ -227,81 +470,90 @@ void discreteDiskFace::buildOct() const
discreteDiskFaceCentroid, discreteDiskFaceInEle);
count = 0;
-
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
MTriangle *t = discrete_triangles[i];
- std::map<MVertex*,SPoint3>::const_iterator it0 =
- coordinates.find(t->getVertex(0));
- std::map<MVertex*,SPoint3>::const_iterator it1 =
- coordinates.find(t->getVertex(1));
- std::map<MVertex*,SPoint3>::const_iterator it2 =
- coordinates.find(t->getVertex(2));
- _ddft[count].p1 = it0->second;
- _ddft[count].p2 = it1->second;
- _ddft[count].p3 = it2->second;
+ discreteDiskFaceTriangle* my_ddft = &_ddft[count];
+ init_ddft(my_ddft); // #fixme destructor
+ for(unsigned int io=0; io<_N; io++){
+ std::map<MVertex*,SPoint3>::const_iterator it =
+ coordinates.find(t->getVertex(io));
+ my_ddft->p[io] = it->second;
+ }
if(geomType() != GEntity::DiscreteDiskSurface){// CAD
if (t->getVertex(0)->onWhat()->dim() == 2){
- reparamMeshEdgeOnFace(t->getVertex(0), t->getVertex(1), _parent,
- _ddft[count].gfp1, _ddft[count].gfp2);
- reparamMeshEdgeOnFace(t->getVertex(0), t->getVertex(2), _parent,
- _ddft[count].gfp1, _ddft[count].gfp3);
+ for(unsigned int io=1; io<_N; io++)
+ reparamMeshEdgeOnFace(t->getVertex(0), t->getVertex(io), _parent,
+ my_ddft->gfp[0], my_ddft->gfp[io]);
}
else if (t->getVertex(1)->onWhat()->dim() == 2){
- reparamMeshEdgeOnFace(t->getVertex(1), t->getVertex(0), _parent,
- _ddft[count].gfp2, _ddft[count].gfp1);
- reparamMeshEdgeOnFace(t->getVertex(1), t->getVertex(2), _parent,
- _ddft[count].gfp2, _ddft[count].gfp3);
+ for(unsigned int io=0; io<_N; io++)
+ if (io != 1)
+ reparamMeshEdgeOnFace(t->getVertex(1), t->getVertex(io), _parent,
+ my_ddft->gfp[1], my_ddft->gfp[io]);
}
else if (t->getVertex(2)->onWhat()->dim() == 2){
- reparamMeshEdgeOnFace(t->getVertex(2), t->getVertex(0), _parent,
- _ddft[count].gfp3, _ddft[count].gfp1);
- reparamMeshEdgeOnFace(t->getVertex(2), t->getVertex(1), _parent,
- _ddft[count].gfp3, _ddft[count].gfp2);
+ for(unsigned int io=0; io<_N; io++)
+ if (io != 2)
+ reparamMeshEdgeOnFace(t->getVertex(2), t->getVertex(io), _parent,
+ my_ddft->gfp[2], my_ddft->gfp[io]);
+ }
+ else {// quid interior vertices ?
+ for(unsigned int io=0; io<_N; io++)
+ reparamMeshVertexOnFace(t->getVertex(io), _parent, my_ddft->gfp[io]);
}
- else {
- reparamMeshVertexOnFace(t->getVertex(0), _parent, _ddft[count].gfp1);
- reparamMeshVertexOnFace(t->getVertex(1), _parent, _ddft[count].gfp2);
- reparamMeshVertexOnFace(t->getVertex(2), _parent, _ddft[count].gfp3);
+ }
+ for(unsigned int io=0; io<_N; io++)
+ my_ddft->v[io] = SPoint3(t->getVertex(io)->x(), t->getVertex(io)->y(),
+ t->getVertex(io)->z());
+
+ my_ddft->gf = _parent;
+ my_ddft->tri = t;
+
+ SVector3 dXdxi (my_ddft->v[1] - my_ddft->v[0]); // constant
+ SVector3 dXdeta(my_ddft->v[2] - my_ddft->v[0]); // constant
+
+ double grads[_N][3];// u,v,w=0,grads
+ for(unsigned int myI=0; myI<t->getNumVertices(); myI++){
+ double U,V,W;
+ t->getNode(myI,U,V,W); // Xi Eta
+ mynodalbasis->df(U,V,0,grads);
+
+ //compute first derivatives for every triangle
+ // (dXi/dU)^-1 = dU/dXi
+ double mat[2][2] = {{0.,0.},{0.,0.}};
+ for(unsigned int io=0; io<_N; io++){
+ for (unsigned int ic=0; ic<2; ic++){
+ mat[0][ic] += my_ddft->p[io].x() * grads[io][ic];
+ mat[1][ic] += my_ddft->p[io].y() * grads[io][ic];
+ }
}
+
+ double inv[2][2];
+ inv2x2(mat,inv);
+
+ SVector3 dXdu(dXdxi * inv[0][0] + dXdeta * inv[1][0]);
+ SVector3 dXdv(dXdxi * inv[0][1] + dXdeta * inv[1][1]);
+ firstElemDerivatives[(MElement*)t].push_back(Pair<SVector3,SVector3>(dXdu,dXdv));
}
- _ddft[count].v1 = SPoint3(t->getVertex(0)->x(), t->getVertex(0)->y(),
- t->getVertex(0)->z());
- _ddft[count].v2 = SPoint3(t->getVertex(1)->x(), t->getVertex(1)->y(),
- t->getVertex(1)->z());
- _ddft[count].v3 = SPoint3(t->getVertex(2)->x(), t->getVertex(2)->y(),
- t->getVertex(2)->z());
- _ddft[count].gf = _parent;
- _ddft[count].tri = t;
-
-
- //compute first derivatives for every triangle , to disappear (erase)
- double mat[2][2] = {{_ddft[count].p2.x() - _ddft[count].p1.x(),
- _ddft[count].p3.x() - _ddft[count].p1.x()},
- {_ddft[count].p2.y() - _ddft[count].p1.y(),
- _ddft[count].p3.y() - _ddft[count].p1.y()}}; // modified higher
-
- double inv[2][2];
- inv2x2(mat, inv);
- SVector3 dXdxi (_ddft[count].v2 - _ddft[count].v1); // constant
- SVector3 dXdeta(_ddft[count].v3 - _ddft[count].v1); // constant
- SVector3 dXdu(dXdxi * inv[0][0] + dXdeta * inv[1][0]); // to be modified for higher order
- SVector3 dXdv(dXdxi * inv[0][1] + dXdeta * inv[1][1]); // to be modified for higher order
- firstElemDerivatives[(MElement*)t] = Pair<SVector3,SVector3>(dXdu,dXdv);
-
- nodes[count][0] = (it0->second.x() + it1->second.x() + it2->second.x())/3.0 ;
- nodes[count][1] = (it0->second.y() + it1->second.y() + it2->second.y())/3.0 ;
- nodes[count][2] = 0.0;
-
- Octree_Insert(&_ddft[count], oct);
+ nodes[count][0] = 0.;
+ nodes[count][1] = 0.;
+ nodes[count][2] = 0.;
+ for(unsigned int io=0; io<_N; io++){
+ nodes[count][0] += my_ddft->p[io].x();
+ nodes[count][1] += my_ddft->p[io].y();
+ }
+ nodes[count][0] /= (double) _N;// #notsure
+ nodes[count][1] /= (double) _N;
+
+ Octree_Insert(my_ddft, oct);
count++;
- }
+ }// end big for over mesh[i]
Octree_Arrange(oct);
- // USELESS, laplacian
//smooth first derivatives at vertices
if(adjv.size() == 0){
- std::vector<MTriangle*> allTri;
+ std::vector<MTriangle*> allTri; //
allTri.insert(allTri.end(), discrete_triangles.begin(), discrete_triangles.end() );
buildVertexToTriangle(allTri, adjv);
}
@@ -311,15 +563,16 @@ void discreteDiskFace::buildOct() const
SVector3 dXdu(0.0), dXdv(0.0);
int nbTri = vTri.size();
for (int j = 0; j < nbTri; j++){ // elements's contribution for a vertex
- dXdu += firstElemDerivatives[vTri[j]].first();
- dXdv += firstElemDerivatives[vTri[j]].second();
+ std::vector<Pair<SVector3,SVector3> > myVec = firstElemDerivatives[vTri[j]];
+ dXdu += myVec[nodeLocalNum(vTri[j],v)].first();
+ dXdv += myVec[nodeLocalNum(vTri[j],v)].second();
}
dXdu *= 1./nbTri;
dXdv *= 1./nbTri;
firstDerivatives[v] = Pair<SVector3, SVector3>(dXdu, dXdv);
}
- kdtree = new ANNkd_tree(nodes, count, 3);
+ kdtree = new ANNkd_tree(nodes, count, _N); // #notsure
}
@@ -347,14 +600,10 @@ bool discreteDiskFace::parametrize() const
for(size_t i = 0; i < discrete_triangles.size(); ++i){
MTriangle *t = discrete_triangles[i];
-
- myAssemblerU.numberVertex(t->getVertex(0), 0, 1);
- myAssemblerU.numberVertex(t->getVertex(1), 0, 1);
- myAssemblerU.numberVertex(t->getVertex(2), 0, 1);
-
- myAssemblerV.numberVertex(t->getVertex(0), 0, 1);
- myAssemblerV.numberVertex(t->getVertex(1), 0, 1);
- myAssemblerV.numberVertex(t->getVertex(2), 0, 1);
+ for(unsigned int j=0; j<t->getNumVertices(); j++){
+ myAssemblerU.numberVertex(t->getVertex(j), 0, 1);
+ myAssemblerV.numberVertex(t->getVertex(j), 0, 1);
+ }
}
@@ -372,10 +621,7 @@ bool discreteDiskFace::parametrize() const
SElement se(discrete_triangles[i]);
mappingU.addToMatrix(myAssemblerU, &se);
mappingV.addToMatrix(myAssemblerV, &se);
- }
-
-
- // lsys_v->printMatlab("testv.m");
+ }
double t2 = Cpu();
@@ -396,7 +642,7 @@ bool discreteDiskFace::parametrize() const
p[1] = value_V;
coordinates[v] = p;
}
- else{
+ else{ // #useful ?
itf->second[0]= value_U;
itf->second[1]= value_V;
}
@@ -414,39 +660,76 @@ void discreteDiskFace::getTriangleUV(const double u,const double v,
double &_u, double &_v)const{
double uv[3] = {u,v,0.};
*mt = (discreteDiskFaceTriangle*) Octree_Search(uv,oct);
- if (!(*mt))return;
-
- double M[2][2],X[2],R[2];
- const SPoint3 p0 = (*mt)->p1;
- const SPoint3 p1 = (*mt)->p2;
- const SPoint3 p2 = (*mt)->p3;
- 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 - p0.x());
- R[1] = (v - p0.y());
- sys2x2(M, R, X);
- _u = X[0];
- _v = X[1];
+ if (!(*mt)) return;
+
+ if (_order == 1){
+ double M[2][2],X[2],R[2];
+ const SPoint3 p0 = (*mt)->p[0];
+ const SPoint3 p1 = (*mt)->p[1];
+ const SPoint3 p2 = (*mt)->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 - p0.x());
+ R[1] = (v - p0.y());
+ sys2x2(M, R, X);
+ _u = X[0];// xi
+ _v = X[1];// eta
+ }
+ else{// newton raphson
+ discreteDiskFaceTriangle *my_ddft = *mt;
+ double Xi[3];
+ my_ddft->tri->xyz2uvw(uv,Xi);
+ _u = Xi[0];
+ _v = Xi[1];
+ }
+}
+void discreteDiskFace::checkOrientationUV(){
+
+ double initial, current; // initial and current orientation
+ discreteDiskFaceTriangle *ct;
+ ct = &_ddft[0];
+ double p1[2] = {ct->p[0].x(), ct->p[0].y()};
+ double p2[2] = {ct->p[1].x(), ct->p[1].y()};
+ double p3[2] = {ct->p[2].x(), ct->p[2].y()};
+ initial = robustPredicates::orient2d(p1, p2, p3);
+ unsigned int i=1;
+ for (; i<discrete_triangles.size(); i++){
+ ct = &_ddft[i];
+ p1[0] = ct->p[0].x(); p1[1] = ct->p[0].y();
+ p2[0] = ct->p[1].x(); p2[1] = ct->p[1].y();
+ p3[0] = ct->p[2].x(); p3[1] = ct->p[2].y();
+ current = robustPredicates::orient2d(p1, p2, p3);
+ if(initial*current < 0.) break;
+ }
+ if(i < discrete_triangles.size())
+ Msg::Error("discreteDiskFace: The parametrization is not one-to-one :-(");
+ else
+ std::cout << "discreteDiskFace: The parametrization is one-to-one :-)" << std::endl;
+
}
// (u;v) |-> < (x;y;z); GFace; (u;v) >
GPoint discreteDiskFace::point(const double par1, const double par2) const
{
-
double U,V;
double par[2] = {par1,par2};
discreteDiskFaceTriangle* dt;
- getTriangleUV(par1,par2,&dt,U,V);
+ getTriangleUV(par1,par2,&dt,U,V);// return Xi,Eta
if (!dt) {
GPoint gp = GPoint(1.e22,1.e22,1.e22,_parent,par);
gp.setNoSuccess();
return gp;
}
- SPoint3 p = dt->v1*(1.-U-V) + dt->v2*U + dt->v3*V;
+ //polynomialBasis myPolynomial = polynomialBasis(dt->tri->getTypeForMSH());
+ double eval[_N];
+ mynodalbasis->f(U,V,0.,eval);//#mark
+ SPoint3 p;
+ for(int io=0; io<_N; io++)
+ p += dt->v[io]*eval[io];
return GPoint(p.x(),p.y(),p.z(),_parent,par);
}
@@ -469,13 +752,13 @@ SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
if (de){
MVertex *v1,*v2;
double xi;
- de->interpolateInGeometry (v,&v1,&v2,xi);
+ de->interpolateInGeometry (v,&v1,&v2,xi); // modify
std::map<MVertex*,SPoint3>::iterator it1 = coordinates.find(v1);
std::map<MVertex*,SPoint3>::iterator it2 = coordinates.find(v2);
if(it1 == coordinates.end()) Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__);
if(it2 == coordinates.end()) Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__);
return SPoint2(it1->second.x(),it1->second.y()) * (1.-xi) +
- SPoint2(it2->second.x(),it2->second.y()) * xi;
+ SPoint2(it2->second.x(),it2->second.y()) * xi; // modify
}
}
Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__);
@@ -506,7 +789,6 @@ double discreteDiskFace::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVec
Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 ¶m) const
{
-
double U,V;
discreteDiskFaceTriangle* ddft;
getTriangleUV(param.x(),param.y(),&ddft,U,V);
@@ -517,34 +799,26 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 ¶m) const
Msg::Warning("discreteDiskFace::firstDer << triangle not found %g %g",param[0],param[1]);
return Pair<SVector3, SVector3>(SVector3(1,0,0), SVector3(0,1,0));
}
+
+ SVector3 dXdU = 0.;//dXdu1*(1.-U-V) + dXdu2*U + dXdu3*V;
+ SVector3 dXdV = 0.;//dXdv1*(1.-U-V) + dXdv2*U + dXdv3*V;
+
+ double eval[_N];
+ mynodalbasis->f(U,V,0.,eval); //#mark
+
+ for (unsigned int io=0; io<_N; io++){
+ std::map<MVertex*, Pair<SVector3,SVector3> >::iterator it =
+ firstDerivatives.find(tri->getVertex(io));
+ if (it == firstDerivatives.end()) Msg::Fatal ("Vertex %d (%d) has no firstDerivatives",
+ tri->getVertex(io)->getNum(),io);
+ SVector3 dXdu = it->second.first();
+ SVector3 dXdv = it->second.second();
+
+ dXdU += dXdu * eval[io];
+ dXdV += dXdv * eval[io];
+ }
- std::map<MVertex*, Pair<SVector3,SVector3> >::iterator it0 =
- firstDerivatives.find(tri->getVertex(0));
- if (it0 == firstDerivatives.end()) Msg::Fatal ("Vertex %d (0) has no firstDerivatives",
- tri->getVertex(0)->getNum());
- std::map<MVertex*, Pair<SVector3,SVector3> >::iterator it1 =
- firstDerivatives.find(tri->getVertex(1));
- if (it1 == firstDerivatives.end()) Msg::Fatal ("Vertex %d (1) has no firstDerivatives",
- tri->getVertex(1)->getNum());
- std::map<MVertex*, Pair<SVector3,SVector3> >::iterator it2 =
- firstDerivatives.find(tri->getVertex(2));
- if (it2 == firstDerivatives.end()) Msg::Fatal ("Vertex %d (2) has no firstDerivatives",
- tri->getVertex(2)->getNum());
-
-
- SVector3 dXdu1 = it0->second.first();
- SVector3 dXdu2 = it1->second.first();
- SVector3 dXdu3 = it2->second.first();
-
- SVector3 dXdv1 = it0->second.second();
- SVector3 dXdv2 = it1->second.second();
- SVector3 dXdv3 = it2->second.second();
-
- SVector3 dXdu = dXdu1*(1.-U-V) + dXdu2*U + dXdu3*V;
- SVector3 dXdv = dXdv1*(1.-U-V) + dXdv2*U + dXdv3*V;
-
-
- return Pair<SVector3, SVector3>(dXdu,dXdv);
+ return Pair<SVector3, SVector3>(dXdU,dXdV);
}
@@ -572,17 +846,64 @@ void discreteDiskFace::buildAllNodes()
void discreteDiskFace::putOnView()
{
+
+
+
+ FILE* view_u = Fopen("param_u.pos","w");
+ FILE* view_v = Fopen("param_v.pos","w");
+ FILE* UVxyz = Fopen("UVxyz.pos","w");
+ if(view_u && view_v){
+ fprintf(view_u,"View \"paramU\"{\n");
+ fprintf(view_v,"View \"paramV\"{\n");
+ fprintf(UVxyz,"View \"Uxyz\"{\n");
+ for (unsigned int i=0; i<discrete_triangles.size(); i++){
+ discreteDiskFaceTriangle* my_ddft = &_ddft[i];
+ if (_order == 1){
+ fprintf(view_u,"ST(");
+ fprintf(view_v,"ST(");
+ fprintf(UVxyz,"ST(");
+ }
+ else{
+ fprintf(view_u,"ST%d(",_order);
+ fprintf(view_v,"ST%d(",_order);
+ fprintf(UVxyz,"ST%d(",_order);
+ }
+ for (unsigned int j=0; j<_N-1; j++){
+ fprintf(view_u,"%g,%g,%g,",my_ddft->v[j].x(),my_ddft->v[j].y(),my_ddft->v[j].z());
+ fprintf(view_v,"%g,%g,%g,",my_ddft->v[j].x(),my_ddft->v[j].y(),my_ddft->v[j].z());
+ fprintf(UVxyz,"%g,%g,%g,",my_ddft->p[j].x(),my_ddft->p[j].y(),0.);
+ }
+ fprintf(view_u,"%g,%g,%g){",my_ddft->v[_N-1].x(),my_ddft->v[_N-1].y(),my_ddft->v[_N-1].z());
+ fprintf(view_v,"%g,%g,%g){",my_ddft->v[_N-1].x(),my_ddft->v[_N-1].y(),my_ddft->v[_N-1].z());
+ fprintf(UVxyz,"%g,%g,%g){",my_ddft->p[_N-1].x(),my_ddft->p[_N-1].y(),0.);
+ for (unsigned int j=0; j<_N-1; j++){
+ fprintf(view_u,"%g,",my_ddft->p[j].x());
+ fprintf(view_v,"%g,",my_ddft->p[j].y());
+ fprintf(UVxyz,"%g,",sqrt(my_ddft->v[j].x()*my_ddft->v[j].x()+my_ddft->v[j].z()*my_ddft->v[j].z()+my_ddft->v[j].y()*my_ddft->v[j].y()));
+ }
+ fprintf(view_u,"%g};\n",my_ddft->p[_N-1].x());
+ fprintf(view_v,"%g};\n",my_ddft->p[_N-1].y());
+ fprintf(UVxyz,"%g};\n",sqrt(my_ddft->v[_N-1].x()*my_ddft->v[_N-1].x()+my_ddft->v[_N-1].z()*my_ddft->v[_N-1].z()+my_ddft->v[_N-1].y()*my_ddft->v[_N-1].y()));
+ }
+ fprintf(view_u,"};\n");
+ fprintf(view_v,"};\n");
+ fprintf(UVxyz,"};\n");
+ fclose(view_u);
+ fclose(view_v);
+ fclose(UVxyz);
+ }
+
+ /*
#ifdef HAVE_POST
std::map<int, std::vector<double> > u;
std::map<int, std::vector<double> > v;
-
for(std::set<MVertex*>::iterator iv = allNodes.begin(); iv!=allNodes.end(); ++iv){
MVertex *p = *iv;
u[p->getNum()].push_back(coordinates[p].x());
v[p->getNum()].push_back(coordinates[p].y());
}
PView* view_u = new PView("u", "NodeData", GModel::current(), u);
- PView* view_v = new PView("u", "NodeData", GModel::current(), v);
+ PView* view_v = new PView("v", "NodeData", GModel::current(), v);
view_u->setChanged(true);
view_v->setChanged(true);
view_u->write("param_u.msh", 5);
@@ -592,6 +913,7 @@ void discreteDiskFace::putOnView()
#else
Msg::Error("discreteDiskFace: cannot export without post module")
#endif
+ */
}
// useful for mesh generators ----------------------------------------
@@ -614,10 +936,8 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
// the point is situated in the half plane defined
// by direction n1 and point p (exclude the one on the
// other side)
-
- SVector3 t1 = ct->v2 - ct->v1;
- SVector3 t2 = ct->v3 - ct->v1;
-
+ SVector3 t1 = ct->v[1] - ct->v[0];
+ SVector3 t2 = ct->v[2] - ct->v[0];
// let us first compute point q to be the intersection
// of the plane of the triangle with the line L = p + t n1
// Compute w = t1 x t2 = (a,b,c) and write a point of the plabe as
@@ -626,15 +946,13 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
// a (p_x + t n1_x) + a (p_y + t n1_y) + c (p_z + t n1_z) -
// (v1_x a + v1_y b + v1_z * c) = 0
// gives t
-
SVector3 w = crossprod(t1,t2);
double t = d;
// if everything is not coplanar ...
if (fabs(dot(n1,w)) > 1.e-12){
- t = dot(ct->v1-p,w)/dot(n1,w);
+ t = dot(ct->v[0]-p,w)/dot(n1,w);
}
SVector3 q = p + n1*t;
-
// consider the line that intersects both planes in its
// parametric form
// X(x,y,z) = q + t m
@@ -650,12 +968,10 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
// we choose the one that corresponds to (s_i - p) . n1 > 0
SVector3 m = crossprod(w,n);
-
const double a = dot(m,m);
const double b = 2*dot(m,q-p);
const double c = dot(q-p,q-p) - d*d;
const double delta = b*b-4*a*c;
-
if (delta >= 0){
const double ta = (-b + sqrt(delta)) / (2.*a);
const double tb = (-b - sqrt(delta)) / (2.*a);
@@ -681,20 +997,20 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
mat[0][0] = dot(t1,t1);
mat[1][1] = dot(t2,t2);
mat[0][1] = mat[1][0] = dot(t1,t2);
- b[0] = dot(s-ct->v1,t1) ;
- b[1] = dot(s-ct->v1,t2) ;
+ b[0] = dot(s-ct->v[0],t1) ;
+ b[1] = dot(s-ct->v[0],t2) ;
sys2x2(mat,b,uv);
// check now if the point is inside the triangle
if (uv[0] >= -1.e-6 && uv[1] >= -1.e-6 &&
1.-uv[0]-uv[1] >= -1.e-6 ) {
SVector3 pp =
- ct->p1 * ( 1.-uv[0]-uv[1] ) +
- ct->p2 *uv[0] +
- ct->p3 *uv[1] ;
+ ct->p[0] * ( 1.-uv[0]-uv[1] ) +
+ ct->p[1] * uv[0] +
+ ct->p[2] * uv[1] ;
return GPoint(s.x(),s.y(),s.z(),this,pp);
}
}
- }
+ }
GPoint pp(0);
pp.setNoSuccess();
Msg::Debug("ARGG no success intersection circle");
@@ -702,4 +1018,6 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
}
+
+
#endif
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index 6aa25db7a7231383a4a4e342e51515c00ccbb075..ee37d5959837b3b8d2042440df3409d68c03c622 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -21,6 +21,7 @@
#include "MElementOctree.h"
#include "meshGFaceOptimize.h"
#include "PView.h"
+#include "robustPredicates.h"
class ANNkd_tree;
class Octree;
@@ -29,9 +30,9 @@ class GRbf;
class discreteDiskFaceTriangle {
public:
- SPoint3 p1, p2, p3; // vertices in (u;v)
- SPoint2 gfp1, gfp2, gfp3; // CAD model
- SPoint3 v1, v2, v3; // vertices in (x;y;z)
+ SPoint3* p; // vertices in (u;v)
+ SPoint2* gfp; // CAD model
+ SPoint3* v;// vertices in (x;y;z)
GFace *gf; // GFace tag
MTriangle *tri; // mesh triangle in (x;y;z)
discreteDiskFaceTriangle() : gf(0), tri(0) {}
@@ -41,14 +42,14 @@ class discreteDiskFaceTriangle {
class discreteDiskFace : public GFace {
GFace *_parent;
- // bool checkOrientation(int iter, bool moveBoundaries=false) const; // mark todo
void buildOct() const;
bool parametrize() const;
void buildAllNodes() ;
void getBoundingEdges();
void putOnView();
+
public:
- discreteDiskFace(GFace *parent, std::vector<MTriangle*> &mesh);
+ discreteDiskFace(GFace *parent, std::vector<MTriangle*> &mesh, int p);// MTriangle -> MTriangle 6
virtual ~discreteDiskFace() {triangles.clear();}
void getTriangleUV(const double u,const double v,discreteDiskFaceTriangle **mt, double &_u, double &_v) const;
GPoint point(double par1, double par2) const;
@@ -63,12 +64,31 @@ class discreteDiskFace : public GFace {
GPoint intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
const SVector3 &p, const double &d,
double uv[2]) const;
+ void checkOrientationUV();
protected:
//------------------------------------------------
// a copy of the mesh that should not be destroyed
- std::vector<MTriangle*> discrete_triangles;
+ std::vector<MTriangle*> discrete_triangles; // MTriangleN AND MTriangle6 are children of MTriangle
std::vector<MVertex*> discrete_vertices;
//------------------------------------------------
+ int nodeLocalNum(MElement* e, MVertex* v) const{
+ for(unsigned int i=0; i<e->getNumVertices(); i++)
+ if (v == e->getVertex(i))
+ return i;
+ };
+ int edgeLocalNum(MElement* e, MEdge ed) const{
+ for(unsigned int i=0; i<e->getNumEdges(); i++)
+ if (ed == e->getEdge(i))
+ return i;
+ };
+ void init_ddft(discreteDiskFaceTriangle* ddft) const{
+ ddft->p = new SPoint3[_N];
+ ddft->gfp = new SPoint2[_N];
+ ddft->v = new SPoint3[_N];
+ };
+ //------------------------------------------------
+ int _order;
+ int _N;
double _totLength;
std::map<double,std::vector<MVertex*> > _loops;
std::vector<MVertex*> _U0; // dirichlet's bc's
@@ -77,13 +97,14 @@ class discreteDiskFace : public GFace {
mutable std::map<SPoint3,SPoint3 > _coordPoints; // ?
mutable v2t_cont adjv; // ? v2t_cont ? ?????
mutable std::map<MVertex*, Pair<SVector3,SVector3> > firstDerivatives;
- mutable std::map<MElement*, Pair<SVector3,SVector3> > firstElemDerivatives;
+ mutable std::map<MElement*, std::vector<Pair<SVector3,SVector3> > > firstElemDerivatives;//dXdU
mutable discreteDiskFaceTriangle *_ddft;
mutable ANNkd_tree *uv_kdtree;
mutable ANNkd_tree *kdtree;
mutable Octree *oct;
mutable std::vector<double> _coords;
-
+ const nodalBasis* mynodalbasis;
+ const bezierBasis* mybezierbasis;
};
#endif
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 1211b07f71f84f23cbb9b55762323197fd6e25a3..fb8388ca6c0ada40c40e6c32ecdc2f042e052233 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -16,7 +16,7 @@
#include "OS.h"
discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
-{
+{
Surface *s = Create_Surface(num, MSH_SURF_DISCRETE);
Tree_Add(model->getGEOInternals()->Surfaces, &s);
meshStatistics.status = GFace::DONE;
@@ -72,79 +72,52 @@ GPoint discreteFace::point(double par1, double par2) const
SPoint2 discreteFace::parFromPoint(const SPoint3 &p, bool onSurface) const
{
- if (getCompound()){
- return getCompound()->parFromPoint(p);
- }
- else{
- Msg::Error("Cannot compute parametric coordinates on discrete face");
- return SPoint2();
- }
+ return SPoint2();
}
SVector3 discreteFace::normal(const SPoint2 ¶m) const
{
- if (getCompound()){
- return getCompound()->normal(param);
- }
- else{
- Msg::Error("Cannot evaluate normal on discrete face");
return SVector3();
- }
}
double discreteFace::curvatureMax(const SPoint2 ¶m) const
{
- if (getCompound()){
- return getCompound()->curvatureMax(param);
- }
- else{
- Msg::Error("Cannot evaluate curvature on discrete face");
+
return false;
- }
}
double discreteFace::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
double *curvMax, double *curvMin) const
{
- if (getCompound()){
- return getCompound()->curvatures(param, dirMax, dirMin, curvMax, curvMin);
- }
- else{
- Msg::Error("Cannot evaluate curvatures and curvature directions on discrete face");
- return false;
- }
+ return false;
}
Pair<SVector3, SVector3> discreteFace::firstDer(const SPoint2 ¶m) const
{
- if (getCompound()){
- return getCompound()->firstDer(param);
- }
- else{
- Msg::Error("Cannot evaluate derivative on discrete face");
- return Pair<SVector3, SVector3>();
- }
+ return Pair<SVector3, SVector3>();
}
void discreteFace::secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
{
+ /*<<<<<<< .mine
+=======
if (getCompound()){
return getCompound()->secondDer(param, dudu, dvdv, dudv);
}
else{
Msg::Error("Cannot evaluate second derivative on discrete face");
+ >>>>>>> .r22873*/
return;
- }
}
// FIXME PAB ----> really create an atlas !!!!!!!!!!!!!!!
void discreteFace::createGeometry()
{
#if defined(HAVE_ANN) && defined(HAVE_SOLVER)
- if (!_atlas.empty())return;
+ if (!_atlas.empty())return;
// parametrization is done here !!!
- discreteDiskFace *df = new discreteDiskFace (this, triangles);
+ discreteDiskFace *df = new discreteDiskFace (this, triangles,2);
df->replaceEdges(l_edges);
_atlas.push_back(df);
#endif