Skip to content
Snippets Groups Projects
Commit d5f0ffd9 authored by PA Beaufort's avatar PA Beaufort
Browse files 

crossField :-)
parent cda67987
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/discreteFace.cpp
+ 341
− 7
View file @ d5f0ffd9
......@@ -13,6 +13,13 @@
#include "OS.h"
#include <stack>
#include <queue>
#include <complex>
// #include <cmath>
#if defined(HAVE_PETSC)
#include "linearSystemPETSc.h"
#endif
#include "MPoint.h"
......@@ -22,6 +29,57 @@ extern "C" {
}
#endif
static inline double getAlpha(MTriangle* tri, int edj){
double alpha;
if (edj==0)
alpha = 0.;
else{
MVertex *v0, *v1, *v2;
v0 = tri->getVertex(0);
v1 = tri->getVertex(1);
v2 = tri->getVertex(2);
SVector3 a(v1->x()-v0->x(),v1->y()-v0->y(),v1->z()-v0->z());
SVector3 b(v2->x()-v0->x(),v2->y()-v0->y(),v2->z()-v0->z());
SVector3 n = crossprod(a,b); n.normalize();
v0 = tri->getEdge(0).getSortedVertex(0);
v1 = tri->getEdge(0).getSortedVertex(1);
SVector3 U(v1->x()-v0->x(),v1->y()-v0->y(),v1->z()-v0->z());
U.normalize();
SVector3 V = crossprod(n,U); V.normalize();
v0 = tri->getEdge(edj).getSortedVertex(0);
v1 = tri->getEdge(edj).getSortedVertex(1);
SVector3 e(v1->x()-v0->x(),v1->y()-v0->y(),v1->z()-v0->z());
e.normalize();
alpha = std::atan2(dot(e,V),dot(e,U));
}
return alpha;
}
static inline void crouzeixRaviart(const std::vector<double> &U,std::vector<double> &F){
F.resize(3);
double xsi[3] = {0.,1.,0.};
double eta[3] = {0.,0.,1.};
for(int i=0; i<3; i++)
F[i] = U[0] * (1.-2.*eta[i]) + U[1] * (2.*(xsi[i]+eta[i])-1.) + U[2] * (1-2.*xsi[i]);
}
discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
{
Surface *s = Create_Surface(num, MSH_SURF_DISCRETE);
......@@ -113,8 +171,12 @@ void discreteFace::createGeometry()
{
checkAndFixOrientation();
#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
int order = 1;
int order = 2;
int nPart = 2;
double eta = 5/(2.*3.14);
if (!_atlas.empty())return;
......@@ -127,6 +189,8 @@ void discreteFace::createGeometry()
std::vector<MElement*> tem(triangles.begin(),triangles.end());
triangulation* init = new triangulation(-1, tem,this);
allEdg2Tri = init->ed2tri;
/*
toSplit.push(init);
if((toSplit.top())->genus()!=0 || (toSplit.top())->aspectRatio() > eta ||
(toSplit.top())->seamPoint){
......@@ -158,11 +222,7 @@ void discreteFace::createGeometry()
updateTopology(toParam);
for(unsigned int i=0; i<toParam.size(); i++){
/*
char name[256];
sprintf(name,"map%d.pos",i);
toParam[i]->print(name,i);
*/
fillHoles(toParam[i]);
//sprintf(name,"mapFilled%d.pos",i);
//toParam[i]->print(name, toParam[i]->idNum);
......@@ -175,7 +235,8 @@ void discreteFace::createGeometry()
df->replaceEdges(toParam[i]->my_GEdges);
_atlas.push_back(df);
}
*/
crossField();
#endif
}
......@@ -708,6 +769,279 @@ void discreteFace::addTriangle(triangulation* trian, MTriangle* t)
#endif
}
void discreteFace::complex_crossField()
{
// COMPLEX linear system
linearSystem<std::complex<double> > * lsys;
#ifdef HAVE_PETSC
lsys = new linearSystemPETSc<std::complex<double> >;
#else
Msg::Fatal("Petsc is required (we do need complex in discreteFace::crossField())");
#endif
std::complex<double> i1(0,1);
dofManager<std::complex<double> > myAssembler(lsys);
std::map<MEdge,int,Less_Edge> ed2key;
for (unsigned int i=0; i<triangles.size(); i++){
MTriangle* tri = triangles[i];
for(int j=0; j<3; j++){
MEdge ed = tri->getEdge(j);
std::vector<int> iTri = allEdg2Tri[ed];
int mini = iTri[0];
if (iTri.size()>1)
mini = iTri[0] < iTri[1] ? iTri[0] : iTri[1];
else{
double alpha = getAlpha(tri,j);
printf("CL: ed[%f,%f]--[%f,%f]\t theta: %f \n",ed.getSortedVertex(0)->x(),ed.getSortedVertex(0)->y(),ed.getSortedVertex(1)->x(),ed.getSortedVertex(1)->y(),alpha*180./M_PI);
myAssembler.fixDof(3*mini+j,0,std::complex<double>(std::exp(-4.*i1*alpha))); // #tocheck
}
int num,s;
triangles[mini]->getEdgeInfo(ed,num,s);
myAssembler.numberDof(3*mini+num,0);
ed2key[ed] = 3*mini+num;
}// end for j
}// end for unsigned int i
double grad[3][3] = {{0.,-2.,0.},{2.,2.,0.},{-2.,0.,0.}};
for (unsigned int i=0; i<triangles.size(); i++){
MTriangle* tri = triangles[i];
double jac[3][3];
double invjac[3][3];
double dJac = tri->getJacobian(0., 0., 0., jac);
inv3x3(jac, invjac);
for(int j=0; j<3; j++){
double alpha_j = getAlpha(tri,j);
std::complex<double> ej(std::exp(4.*i1*alpha_j));
Dof R(ed2key[tri->getEdge(j)],0);
for(int k=0; k<3; k++){
double alpha_k = getAlpha(tri,k);
std::complex<double> ek(std::exp(-4.*i1*alpha_k));
std::complex<double> K_jk = 0.;
for (int l=0; l<3; l++)
for(int jj=0; jj<3; jj++)
for(int kk=0; kk<3; kk++)
K_jk += grad[j][jj] * invjac[l][jj] * invjac[l][kk] * grad[k][kk];
K_jk *= ej*ek*dJac/2.;
Dof C(ed2key[tri->getEdge(k)],0);
myAssembler.assemble(R,C,K_jk);
}// end for k
}// end for j
}// end for unsigned int
lsys->systemSolve();
FILE* myfile = Fopen("crossField.pos","w");
fprintf(myfile,"View \"cross\"{\n");
for(unsigned int i=0; i<triangles.size(); i++){
fprintf(myfile,"VT(");
MTriangle* tri = triangles[i];
MEdge ed = tri->getEdge(0);
MVertex *v0 = ed.getSortedVertex(0), *v1=ed.getSortedVertex(1);
SVector3 e(v1->x()-v0->x(),v1->y()-v0->y(),v1->z()-v0->z()); e.normalize();
MEdge ed1 = tri->getEdge(1); MVertex *v10 = ed1.getVertex(0), *v11=ed1.getVertex(1);
SVector3 e1(v11->x()-v10->x(),v11->y()-v10->y(),v11->z()-v10->z());
SVector3 n = crossprod(e,e1); n.normalize();
SVector3 d = crossprod(e,n); d.normalize();
std::vector<double> U; U.resize(3);
std::vector<double> V; V.resize(3);
for(int j=0; j<3; j++){
fprintf(myfile,"%f,%f,%f",tri->getVertex(j)->x(),tri->getVertex(j)->y(),tri->getVertex(j)->z());
if (j<2) fprintf(myfile,",");
MEdge ed = tri->getEdge(j);
double alpha = getAlpha(tri,j);
std::complex<double> cross, Cross(std::exp(4.*i1*alpha));
myAssembler.getDofValue(ed2key[ed],0,cross); // conjugate dof in the local edge basis
Cross *= std::conj(cross); // dof of the local tri basis
U[j] = std::real(Cross);
V[j] = std::imag(Cross);
printf("sol: ed%d %f (tri) ~ %f (ed) \n",j,std::atan2(V[j],U[j])*180./(4.*M_PI),(std::atan2(-std::imag(cross),std::real(cross)))*180./(4.*M_PI));
}
fprintf(myfile,")");
std::vector<double> Fu, Fv;
crouzeixRaviart(U,Fu);
crouzeixRaviart(V,Fv);
fprintf(myfile,"{");
for(int j=0; j<3; j++){
double u = Fu[j], v = Fv[j];
printf("sol: ve%d (%f %f)\n",j,u,v);
fprintf(myfile,"%f,%f,%f",u*e.x()+v*d.x(),u*e.y()+v*d.y(),u*e.z()+v*d.z());
if (j<2) fprintf(myfile,",");
}
fprintf(myfile,"};\n");
}
fprintf(myfile,"};");
fclose(myfile);
}
void discreteFace::crossField()
{
// linear system
linearSystem<double> * lsys;
#ifdef HAVE_PETSC
lsys = new linearSystemPETSc<double>;
#else
Msg::Fatal("Petsc is required ");
#endif
dofManager<double> myAssembler(lsys);
std::map<MEdge,int,Less_Edge> ed2key;
for (unsigned int i=0; i<triangles.size(); i++){
MTriangle* tri = triangles[i];
for(int j=0; j<3; j++){
MEdge ed = tri->getEdge(j);
std::vector<int> iTri = allEdg2Tri[ed];
int mini = iTri[0];
if (iTri.size()>1)
mini = iTri[0] < iTri[1] ? iTri[0] : iTri[1];
else{
double alpha = getAlpha(tri,j);
printf("CL: ed[%f,%f]--[%f,%f]\t Theta: %f \n",ed.getSortedVertex(0)->x(),ed.getSortedVertex(0)->y(),ed.getSortedVertex(1)->x(),ed.getSortedVertex(1)->y(),alpha*180./M_PI);
myAssembler.fixDof(3*mini+j,0,1.); // setting theta and not Theta
myAssembler.fixDof(3*mini+j,1,0.);
}
int num,s;
triangles[mini]->getEdgeInfo(ed,num,s);
myAssembler.numberDof(3*mini+num,0); // u, not U
myAssembler.numberDof(3*mini+num,1); // v, not V
ed2key[ed] = 3*mini+num;
}// end for j
}// end for unsigned int i
double grad[3][3] = {{0.,-2.,0.},{2.,2.,0.},{-2.,0.,0.}};
for (unsigned int i=0; i<triangles.size(); i++){
MTriangle* tri = triangles[i];
double jac[3][3];
double invjac[3][3];
double dJac = tri->getJacobian(0., 0., 0., jac);
inv3x3(jac, invjac);
for(int j=0; j<3; j++){
double alpha_j = getAlpha(tri,j);
Dof Ru(ed2key[tri->getEdge(j)],0);
Dof Rv(ed2key[tri->getEdge(j)],1);
for(int k=0; k<3; k++){
double alpha_k = getAlpha(tri,k);
Dof Cu(ed2key[tri->getEdge(k)],0);
Dof Cv(ed2key[tri->getEdge(k)],1);
double K_jk = 0.;
for (int l=0; l<3; l++)
for(int jj=0; jj<3; jj++)
for(int kk=0; kk<3; kk++)
K_jk += grad[j][jj] * invjac[l][jj] * invjac[l][kk] * grad[k][kk];
K_jk *= dJac/2.;
printf("%f\t",K_jk);
//printf("%f \t %f \n %f \t %f \n",cos(4.*(alpha_j-alpha_k))*K_jk,sin(4.*(alpha_j-alpha_k))*K_jk,sin(4.*(alpha_k-alpha_j))*K_jk,cos(4.*(alpha_j-alpha_k))*K_jk);
myAssembler.assemble(Ru,Cu,cos(4.*(alpha_j-alpha_k))*K_jk);
myAssembler.assemble(Ru,Cv,sin(4.*(alpha_j-alpha_k))*K_jk);
myAssembler.assemble(Rv,Cu,sin(4.*(alpha_k-alpha_j))*K_jk);
myAssembler.assemble(Rv,Cv,cos(4.*(alpha_j-alpha_k))*K_jk);
}// end for k
printf("\n");
}// end for j
printf("\n");
}// end for unsigned int
lsys->systemSolve();
FILE* myfile = Fopen("crossField.pos","w");
fprintf(myfile,"View \"cross\"{\n");
for(unsigned int i=0; i<triangles.size(); i++){
fprintf(myfile,"VT(");
MTriangle* tri = triangles[i];
MEdge ed = tri->getEdge(0);
MVertex *v0, *v1, *v2;
v0 = tri->getVertex(0);
v1 = tri->getVertex(1);
v2 = tri->getVertex(2);
SVector3 a(v1->x()-v0->x(),v1->y()-v0->y(),v1->z()-v0->z());
SVector3 b(v2->x()-v0->x(),v2->y()-v0->y(),v2->z()-v0->z());
SVector3 n = crossprod(a,b); n.normalize();
v0 = ed.getSortedVertex(0);
v1 = ed.getSortedVertex(1);
SVector3 e(v1->x()-v0->x(),v1->y()-v0->y(),v1->z()-v0->z());
e.normalize();
printf("e=(%f %f)\n",e.x(),e.y());
SVector3 d = crossprod(n,e); d.normalize();
printf("d=(%f %f)\n",d.x(),d.y());
std::vector<double> U; U.resize(3);
std::vector<double> V; V.resize(3);
for(int j=0; j<3; j++){
fprintf(myfile,"%f,%f,%f",tri->getVertex(j)->x(),tri->getVertex(j)->y(),tri->getVertex(j)->z());
if (j<2) fprintf(myfile,",");
MEdge ed = tri->getEdge(j);
double u, v;// edge basis
myAssembler.getDofValue(ed2key[ed],0,u);// u
myAssembler.getDofValue(ed2key[ed],1,v);// v
double alpha = getAlpha(tri,j);// triangle basis
U[j] = cos(4.*alpha)*u - sin(4.*alpha)*v;// U, not u
V[j] = sin(4.*alpha)*u + cos(4.*alpha)*v;// V, not v
printf("sol: ed%d[%f,%f]--[%f,%f] \t Theta = %f (tri) ~ theta = %f (ed) ~ alpha = %f ~u=%f, v=%f VS U=%f, V=%f \n",j,ed.getSortedVertex(0)->x(),ed.getSortedVertex(0)->y(),ed.getSortedVertex(1)->x(),ed.getSortedVertex(1)->y(),std::atan2(V[j],U[j])*180./(4.*M_PI),std::atan2(v,u)*180./(4.*M_PI),alpha*180./M_PI,u,v,U[j],V[j]);
}
fprintf(myfile,")");
std::vector<double> Fu, Fv;
crouzeixRaviart(U,Fu);
crouzeixRaviart(V,Fv);
fprintf(myfile,"{");
for(int j=0; j<3; j++){
double u = Fu[j], v = Fv[j]; double theta = std::atan2(v,u)/4.;
printf("theta_%d = %f\n",j,theta*180./M_PI);
fprintf(myfile,"%f,%f,%f",cos(theta)*e.x()-sin(theta)*e.y(),sin(theta)*e.x()+cos(theta)*e.y(),e.z());
if (j<2) fprintf(myfile,",");
}
fprintf(myfile,"};\n");
}
fprintf(myfile,"};");
fclose(myfile);
}
void discreteFace::writeGEO(FILE *fp)
{
fprintf(fp, "Discrete Face(%d) = {",tag());
......
Geo/discreteFace.h
+ 3
− 3
View file @ d5f0ffd9
......@@ -31,7 +31,6 @@ class discreteFace : public GFace {
discreteFace(GModel *model, int num);
virtual ~discreteFace() {}
void checkAndFixOrientation();
//void checkConnectivity(std::vector<std::vector<MElement*> >&); undefined
void setupDiscreteVertex(GVertex*,MVertex*,std::set<MVertex*>*);
void setupDiscreteEdge(discreteEdge*,std::vector<MLine*>,std::set<MVertex*>*);
void splitDiscreteEdge(GEdge*,GVertex*,discreteEdge*[2]);
......@@ -39,6 +38,8 @@ class discreteFace : public GFace {
void split(triangulation*,std::vector<triangulation*>&,int);
void fillHoles(triangulation*);
void addTriangle(triangulation*,MTriangle*);
void complex_crossField();
void crossField();
GPoint point(double par1, double par2) const;
SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
SVector3 normal(const SPoint2 &param) const;
......@@ -56,10 +57,9 @@ class discreteFace : public GFace {
void createGeometry();
void gatherMeshes();
virtual void mesh (bool verbose);
//void printAtlasMesh(std::vector<MElement*>, int);
//void printAtlasMesh(discreteDiskFace*, int);
std::vector<discreteDiskFace*> _atlas;
std::vector<GFace*> _CAD;
std::map<MEdge,std::vector<int>,Less_Edge> allEdg2Tri;
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment