Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

source: issm/trunk/src/c/classes/bamg/Geometry.cpp@ 13395

Last change on this file since 13395 was 13395, checked in by Mathieu Morlighem, 12 years ago
merged trunk-jpl and trunk for revision 13393
File size: 30.2 KB

Rev	Line
[3913]	1	#include <cstdio>
	2	#include <cstring>
	3	#include <cmath>
	4	#include <ctime>
	5
	6	#include "../objects.h"
	7	#include "../../include/include.h"
	8	#include "../../shared/Exceptions/exceptions.h"
	9
	10	namespace bamg {
	11
	12	static const Direction NoDirOfSearch=Direction();
	13
	14	/Constructors/Destructors/
[12365]	15	/FUNCTION Geometry::Geometry(){{{/
[5120]	16	Geometry::Geometry(){
	17	Init();
	18	}
[5130]	19	/}}}/
[12365]	20	/FUNCTION Geometry::Geometry(BamgGeom bamggeom, BamgOpts* bamgopts){{{*/
[5120]	21	Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
	22	Init();
	23	ReadGeometry(bamggeom,bamgopts);
[5397]	24	PostRead();
[5120]	25	}
[5130]	26	/}}}/
[12365]	27	/FUNCTION Geometry::Geometry(const Geometry & Gh) (COPY operator){{{/
[3913]	28	Geometry::Geometry(const Geometry & Gh) {
	29	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Geometry)/
	30
	31	long i;
	32	*this = Gh;
	33	NbRef =0;
	34	quadtree=0;
[5573]	35	vertices = nbv ? new GeomVertex[nbv] : NULL;
	36	edges = nbe ? new GeomEdge[nbe]:NULL;
[5148]	37	curves= nbcurves ? new Curve[nbcurves]:NULL;
[5573]	38	subdomains = nbsubdomains ? new GeomSubDomain[nbsubdomains]:NULL;
[3913]	39	for (i=0;i<nbe;i++)
	40	edges[i].Set(Gh.edges[i],Gh,*this);
[5148]	41	for (i=0;i<nbcurves;i++)
[3913]	42	curves[i].Set(Gh.curves[i],Gh,*this);
[5148]	43	for (i=0;i<nbsubdomains;i++)
[3913]	44	subdomains[i].Set(Gh.subdomains[i],Gh,*this);
	45	}
[12365]	46	/}}}/
	47	/FUNCTION Geometry::~Geometry(){{{/
[3913]	48	Geometry::~Geometry() {
	49	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/~Geometry)/
[12509]	50	if(NbRef>0){ _printString_("Trying to delete geometry and NbRef>0, probably due to an error"); return;}
[5208]	51	if(vertices) delete [] vertices; vertices=0;
	52	if(edges) delete [] edges; edges=0;
	53	if(quadtree) delete quadtree; quadtree=0;
	54	if(curves) delete [] curves; curves=0;nbcurves=0;
[3913]	55	if(subdomains) delete [] subdomains;subdomains=0;
[5120]	56	Init();
[3913]	57	}
[12365]	58	/}}}/
[3913]	59
	60	/IO/
[12365]	61	/FUNCTION Geometry::ReadGeometry{{{/
[3913]	62	void Geometry::ReadGeometry(BamgGeom* bamggeom,BamgOpts* bamgopts){
	63
	64	int verbose;
[5120]	65	nbv=0;
	66	nbe=0;
[5148]	67	nbcurves=0;
[3913]	68
	69	double Hmin = HUGE_VAL;// the infinie value
[5143]	70	int i,j,k,n,i0,i1,i2,i3;
[3913]	71
[5120]	72	/initialize some variables/
[5143]	73	verbose= bamgopts->verbose;
	74	nbv = bamggeom->VerticesSize[0];
	75	nbe = bamggeom->EdgesSize[0];
[3913]	76
	77	//some checks
[13036]	78	if (bamggeom->Vertices==NULL) _error_("the domain provided does not contain any vertex");
	79	if (bamggeom->Edges==NULL) _error_("the domain provided does not contain any edge");
[3913]	80
	81	//Vertices
	82	if (bamggeom->Vertices){
[12509]	83	if(verbose>5) _printLine_(" processing Vertices");
[13036]	84	if (bamggeom->VerticesSize[1]!=3) _error_("Vertices should have 3 columns");
[5573]	85	vertices = new GeomVertex[nbv];
[3913]	86	for (i=0;i<nbv;i++) {
	87	vertices[i].r.x=(double)bamggeom->Vertices[i*3+0];
	88	vertices[i].r.y=(double)bamggeom->Vertices[i*3+1];
	89	vertices[i].ReferenceNumber=(long)bamggeom->Vertices[i*3+2];
	90	vertices[i].DirOfSearch=NoDirOfSearch;
	91	vertices[i].color =0;
[5124]	92	vertices[i].type=0;
[3913]	93	}
[5120]	94	/find domain extrema (pmin,pmax) that will define the square box used for by the quadtree/
[3913]	95	pmin = vertices[0].r;
	96	pmax = vertices[0].r;
	97	for (i=0;i<nbv;i++) {
	98	pmin.x = Min(pmin.x,vertices[i].r.x);
	99	pmin.y = Min(pmin.y,vertices[i].r.y);
	100	pmax.x = Max(pmax.x,vertices[i].r.x);
	101	pmax.y = Max(pmax.y,vertices[i].r.y);
	102	}
[5120]	103	/Offset pmin and pmax to avoid round-off errors/
	104	R2 offset = (pmax-pmin)*0.05;
	105	pmin -= offset;
	106	pmax += offset;
	107	/*coefIcoor is the coefficient used for integer coordinates:
	108	* (x-pmin.x)
	109	* Icoor x = (2^30 -1) ------------
	110	* D
	111	* where D is the longest side of the domain (direction x or y)
	112	* so that (x-pmin.x)/D is in ]0 1[
[5177]	113	*
	114	* coefIcoor = (2^30 -1)/D
[5120]	115	*/
	116	coefIcoor=(MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
[13036]	117	if(coefIcoor<=0) _error_("coefIcoor should be positive");
[3913]	118	}
	119	else{
[13036]	120	_error_("No BamgVertex provided");
[3913]	121	}
	122
	123	//Edges
	124	if (bamggeom->Edges){
[5143]	125	R2 zerovector(0,0);
	126	double* verticeslength=NULL;
[3913]	127
[12509]	128	if(verbose>5) _printLine_(" processing Edges");
[13036]	129	if (bamggeom->EdgesSize[1]!=3) _error_("Edges should have 3 columns");
[5573]	130	edges = new GeomEdge[nbe];
[3913]	131
[5143]	132	//initialize verticeslength (sum of the lengths of the edges holding vertex)
	133	verticeslength = new double[nbv];
	134	for(i=0;i<nbv;i++) verticeslength[i]=0;
[3913]	135
[5143]	136	/Loop over the edges/
[3913]	137	for (i=0;i<nbe;i++){
	138
	139	i1=(int)bamggeom->Edges[i*3+0]-1; //-1 for C indexing
	140	i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing
[5143]	141	edges[i].v[0]= vertices + i1; //pointer toward vertex i1 (=&vertices[i1])
	142	edges[i].v[1]= vertices + i2; //pointer toward vertex i2
[5148]	143	edges[i].ReferenceNumber=(long)bamggeom->Edges[i*3+2];
[3913]	144
	145	//get length of edge
[5143]	146	R2 x12=vertices[i2].r-vertices[i1].r;
[3913]	147	double l12=sqrt((x12,x12));
	148	Hmin=Min(Hmin,l12);
	149
	150	//initialize other fields
[5143]	151	edges[i].tg[0]=zerovector;
	152	edges[i].tg[1]=zerovector;
[5401]	153	edges[i].AdjVertexIndex[0] = edges[i].AdjVertexIndex[1] = -1;
[5143]	154	edges[i].Adj[0] = edges[i].Adj[1] = NULL;
[5146]	155	edges[i].type = 0;
[3913]	156
	157	//Cracked?
[5148]	158	if (edges[i].ReferenceNumber!=1) edges[i].SetCracked();
[3913]	159
	160	//prepare metric
	161	vertices[i1].color++;
	162	vertices[i2].color++;
[5143]	163	verticeslength[i1] += l12;
	164	verticeslength[i2] += l12;
[3913]	165	}
	166
[5397]	167	// definition the default of the given mesh size
[5143]	168	for (i=0;i<nbv;i++) {
	169	if (vertices[i].color > 0)
	170	vertices[i].m=Metric(verticeslength[i] /(double) vertices[i].color);
	171	else
	172	vertices[i].m=Metric(Hmin);
	173	}
	174	delete [] verticeslength;
[3913]	175
	176	}
	177	else{
[13036]	178	_error_("No edges provided");
[3913]	179	}
	180
	181	//hVertices
[5187]	182	if(bamgopts->hVertices && bamgopts->hVerticesSize[0]==nbv){
[12509]	183	if(verbose>5) _printLine_(" processing hVertices");
[3913]	184	for (i=0;i< nbv;i++){
[13255]	185	if (!xIsNan<IssmPDouble>(bamgopts->hVertices[i])){
[5187]	186	vertices[i].m=Metric((double)bamgopts->hVertices[i]);
[3913]	187	}
	188	}
	189	}
	190
	191	//MetricVertices
[5187]	192	if(bamgopts->metric && bamgopts->metric[0]==nbv){
[12509]	193	if(verbose>5) _printLine_(" processing MetricVertices");
[3913]	194	for (i=0;i< nbv;i++) {
[5187]	195	vertices[i].m = Metric((double)bamgopts->metric[i3+0],(double)bamgopts->metric[i3+1],(double)bamgopts->metric[i*3+2]);
[3913]	196	}
	197	}
	198
[5091]	199	//MaxCornerAngle
	200	if (bamgopts->MaxCornerAngle){
[12509]	201	if(verbose>5) _printLine_(" processing MaxCornerAngle");
[5091]	202	MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180;
[3913]	203	}
	204
	205	//TangentAtEdges
	206	if (bamggeom->TangentAtEdges){
[12509]	207	if(verbose>5) _printString_(" processing TangentAtEdges");
[13036]	208	if (bamggeom->TangentAtEdgesSize[1]!=4) _error_("TangentAtEdges should have 4 columns");
[3913]	209	int n,i,j,k;
	210	R2 tg;
	211
	212	n=bamggeom->TangentAtEdgesSize[0];
	213	for (k=0;k<n;k++) {
	214	i=(int)bamggeom->TangentAtEdges[k*4+0]-1; //for C indexing
	215	j=(int)bamggeom->TangentAtEdges[k*4+1]-1; //for C indexing
	216	tg.x=bamggeom->TangentAtEdges[k*4+2];
	217	tg.y=bamggeom->TangentAtEdges[k*4+3];
[13036]	218	if (i<0 \|\| i>=nbe) _error_("TangentAtEdges first index exceeds matrix dimension");
	219	if (j!=0 && j!=1) _error_("TangentAtEdges second index should be 1 or 2 only");
[3913]	220	edges[i].tg[j] = tg;
	221	}
	222	}
	223
	224	//Corners
	225	if(bamggeom->Corners){
[12509]	226	if(verbose>5) _printString_(" processing Corners");
[13036]	227	if (bamggeom->CornersSize[1]!=1) _error_("Corners should have 1 column");
[3913]	228	n=bamggeom->CornersSize[0];
	229	for (i=0;i<n;i++) {
	230	j=(int)bamggeom->Corners[i]-1; //for C indexing
[13036]	231	if (j>nbv-1 \|\| j<0) _error_("Bad corner definition: should in [0 " << nbv << "]");
[5143]	232	/Required => at the same time SetRequired and SetCorner/
[3913]	233	vertices[j].SetCorner();
[5143]	234	vertices[j].SetRequired();
	235	}
[3913]	236	}
	237
	238	//RequiredVertices
	239	if(bamggeom->RequiredVertices){
[12509]	240	if(verbose>5) _printLine_(" processing RequiredVertices");
[13036]	241	if (bamggeom->RequiredVerticesSize[1]!=1) _error_("RequiredVertices should have 1 column");
[3913]	242	n=bamggeom->RequiredVerticesSize[0];
	243	for (i=0;i<n;i++) {
	244	j=(int)bamggeom->RequiredVertices[i]-1; //for C indexing
[13036]	245	if (j>nbv-1 \|\| j<0) _error_("Bad RequiredVerticess definition: should in [0 " << nbv << "]");
[5143]	246	vertices[j].SetRequired();
	247	}
[3913]	248	}
	249
	250	//RequiredEdges
	251	if(bamggeom->RequiredEdges){
[12509]	252	if(verbose>5) _printLine_(" processing RequiredEdges");
[13036]	253	if (bamggeom->RequiredEdgesSize[1]!=1) _error_("RequiredEdges should have 1 column");
[3913]	254	n=bamggeom->RequiredEdgesSize[0];
	255	for (i=0;i<n;i++) {
	256	j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing
[13036]	257	if (j>nbe-1 \|\| j<0) _error_("Bad RequiredEdges definition: should in [0 " << nbe << "]");
[3913]	258	edges[j].SetRequired();
	259	}
	260	}
	261
	262	//SubDomain
	263	if(bamggeom->SubDomains){
[12509]	264	if(verbose>5) _printLine_(" processing SubDomains");
[13036]	265	if (bamggeom->SubDomainsSize[1]!=4) _error_("SubDomains should have 4 columns");
[5148]	266	nbsubdomains=bamggeom->SubDomainsSize[0];
[5573]	267	subdomains = new GeomSubDomain[nbsubdomains];
[5148]	268	for (i=0;i<nbsubdomains;i++){
[3913]	269	i0=(int)bamggeom->SubDomains[i*4+0];
	270	i1=(int)bamggeom->SubDomains[i*4+1];
	271	i2=(int)bamggeom->SubDomains[i*4+2];
	272	i3=(int)bamggeom->SubDomains[i*4+3];
[13036]	273	if (i0!=2) _error_("Bad Subdomain definition: first number should be 2 (for Edges)");
	274	if (i1>nbe \|\| i1<=0) _error_("Bad Subdomain definition: second number should in [1 " << nbe << "] (edge number)");
[3913]	275	subdomains[i].edge=edges + (i1-1);
[5148]	276	subdomains[i].direction = (int) i2;
	277	subdomains[i].ReferenceNumber = i3;
[3913]	278	}
	279	}
	280	}
[12365]	281	/}}}/
	282	/FUNCTION Geometry::WriteGeometry{{{/
[3913]	283	void Geometry::WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
	284
	285	int verbose;
	286	int nbreq=0;
	287	int nbreqv=0;
	288	int nbtan=0;
	289	int nbcracked=0;
	290	int i,count;
	291
	292	/Get options/
	293	verbose=bamgopts->verbose;
	294
	295	/Vertices/
[12509]	296	if(verbose>5) _printLine_(" writing Vertices");
[3913]	297	bamggeom->VerticesSize[0]=nbv;
	298	bamggeom->VerticesSize[1]=3;
	299	if (nbv){
[12456]	300	bamggeom->Vertices=xNew<double>(3*nbv);
[3913]	301	for (i=0;i<nbv;i++){
	302	bamggeom->Vertices[i*3+0]=vertices[i].r.x;
	303	bamggeom->Vertices[i*3+1]=vertices[i].r.y;
[5148]	304	bamggeom->Vertices[i*3+2]=vertices[i].GetReferenceNumber();
[3913]	305
	306	//update counters
	307	if (vertices[i].Required()) nbreqv++;
	308	}
	309	}
	310
	311	/Edges/
[12509]	312	if(verbose>5) _printLine_(" writing Edges");
[3913]	313	bamggeom->EdgesSize[0]=nbe;
	314	bamggeom->EdgesSize[1]=3;
	315	if (nbe){
[12456]	316	bamggeom->Edges=xNew<double>(3*nbe);
[3913]	317	for (i=0;i<nbe;i++){
[5149]	318	bamggeom->Edges[i*3+0]=GetId(edges[i][0])+1; //back to Matlab indexing
	319	bamggeom->Edges[i*3+1]=GetId(edges[i][1])+1; //back to Matlab indexing
[5148]	320	bamggeom->Edges[i*3+2]=(double)edges[i].ReferenceNumber;
[3913]	321
	322	//update counters
	323	if (edges[i].Required()) nbreq++;
	324	if (edges[i].TgA() && edges[i][0].Corner()) nbtan++;
	325	if (edges[i].TgB() && edges[i][1].Corner()) nbtan++;
	326	}
	327	}
	328
	329	/RequiredEdges/
[12509]	330	if(verbose>5) _printLine_(" writing " << nbreq << " RequiredEdges");
[3913]	331	bamggeom->RequiredEdgesSize[0]=nbreq;
	332	bamggeom->RequiredEdgesSize[1]=1;
	333	if (nbreq){
[12456]	334	bamggeom->RequiredEdges=xNew<double>(1*nbreq);
[3913]	335	count=0;
	336	for (i=0;i<nbe;i++){
	337	if (edges[i].Required()){
	338	bamggeom->RequiredEdges[count]=i+1; //back to Matlab indexing
	339	count=count+1;
	340	}
	341	}
	342	}
	343
	344	//No corners
	345
	346	/RequiredVertices/
[12509]	347	if(verbose>5) _printLine_(" writing " << nbreqv << " RequiredVertices");
[3913]	348	bamggeom->RequiredVerticesSize[0]=nbreqv;
	349	bamggeom->RequiredVerticesSize[1]=1;
	350	if (nbreqv){
[12456]	351	bamggeom->RequiredVertices=xNew<double>(1*nbreqv);
[3913]	352	count=0;
	353	for (i=0;i<nbv;i++){
	354	if (vertices[i].Required()){
	355	bamggeom->RequiredVertices[count]=i+1; //back to Matlab indexing
	356	count=count+1;
	357	}
	358	}
	359	}
	360
	361	/SubDomains/
[12509]	362	if(verbose>5) _printLine_(" writing SubDomains");
[5148]	363	bamggeom->SubDomainsSize[0]=nbsubdomains;
[3913]	364	bamggeom->SubDomainsSize[1]=4;
[5148]	365	if (nbsubdomains){
[12456]	366	bamggeom->SubDomains=xNew<double>(4*nbsubdomains);
[5148]	367	for (i=0;i<nbsubdomains;i++){
[3913]	368	bamggeom->SubDomains[4*i+0]=2;
[5149]	369	bamggeom->SubDomains[4*i+1]=GetId(subdomains[i].edge)+1; //back to Matlab indexing
[5148]	370	bamggeom->SubDomains[4*i+2]=subdomains[i].direction;
	371	bamggeom->SubDomains[4*i+3]=subdomains[i].ReferenceNumber;
[3913]	372	}
	373	}
	374
	375	/TangentAtEdges/
[12509]	376	if(verbose>5) _printLine_(" writing TangentAtEdges");
[3913]	377	bamggeom->TangentAtEdgesSize[0]=nbtan;
	378	bamggeom->TangentAtEdgesSize[1]=4;
	379	if (nbtan){
[12456]	380	bamggeom->TangentAtEdges=xNew<double>(4*nbtan);
[3913]	381	count=0;
	382	for (i=0;i<nbe;i++){
	383	if (edges[i].TgA() && edges[i][0].Corner()){
	384	bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
	385	bamggeom->TangentAtEdges[4*i+1]=1;
	386	bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[0].x;
	387	bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[0].y;
	388	}
	389	if (edges[i].TgB() && edges[i][1].Corner()){
	390	bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
	391	bamggeom->TangentAtEdges[4*i+1]=2;
	392	bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[1].x;
[5016]	393	bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[1].y;
[3913]	394	}
	395	count=count+1;
	396	}
	397	}
	398	}
[12365]	399	/}}}/
[3913]	400
	401	/Methods/
[12365]	402	/FUNCTION Geometry::Echo {{{/
[5400]	403	void Geometry::Echo(void){
	404
[12509]	405	_printLine_("Geometry:");
	406	_printLine_(" nbv (number of vertices) : " << nbv);
	407	_printLine_(" nbe (number of edges) : " << nbe);
	408	_printLine_(" nbsubdomains: " << nbsubdomains);
	409	_printLine_(" nbcurves: " << nbcurves);
	410	_printLine_(" vertices: " << vertices);
	411	_printLine_(" edges: " << edges);
	412	_printLine_(" quadtree: " << quadtree);
	413	_printLine_(" subdomains: " << subdomains);
	414	_printLine_(" curves: " << curves);
	415	_printLine_(" pmin (x,y): (" << pmin.x << " " << pmin.y << ")");
	416	_printLine_(" pmax (x,y): (" << pmax.x << " " << pmax.y << ")");
	417	_printLine_(" coefIcoor: " << coefIcoor);
	418	_printLine_(" MaxCornerAngle: " << MaxCornerAngle);
[5400]	419
	420	return;
	421	}
	422	/}}}/
[12365]	423	/FUNCTION Geometry::Init{{{/
[5400]	424	void Geometry::Init(void){
	425	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)/
	426
	427	NbRef=0;
	428	nbv=0;
	429	nbe=0;
	430	quadtree=NULL;
	431	curves=NULL;
	432	edges=NULL;
	433	vertices=NULL;
	434	nbsubdomains=0;
	435	nbcurves=0;
	436	subdomains=NULL;
	437	MaxCornerAngle = 10*Pi/180; //default is 10 degres
	438	}
[12365]	439	/}}}/
	440	/FUNCTION Geometry::MinimalHmin{{{/
[5400]	441	double Geometry::MinimalHmin() {
	442	/* coeffIcoor = (2^30-1)/D
	443	* We cannot go beyond hmin = D/2^30 because of the quadtree
	444	* hmin is therefore approximately 2/coeffIcoor */
	445	return 2.0/coefIcoor;
	446	}/}}}/
[12365]	447	/FUNCTION Geometry::MaximalHmax{{{/
[5400]	448	double Geometry::MaximalHmax() {
	449	return Max(pmax.x-pmin.x,pmax.y-pmin.y);
	450	}/}}}/
[12365]	451	/FUNCTION Geometry::GetId(const GeomVertex &t){{{/
[5573]	452	long Geometry::GetId(const GeomVertex & t) const {
[5400]	453	return &t - vertices;
	454	}/}}}/
[12365]	455	/FUNCTION Geometry::GetId(const GeomVertex t){{{*/
[5573]	456	long Geometry::GetId(const GeomVertex * t) const {
[5400]	457	return t - vertices;
	458	}/}}}/
[12365]	459	/FUNCTION Geometry::GetId(const GeomEdge & t){{{/
[5573]	460	long Geometry::GetId(const GeomEdge & t) const {
[5400]	461	return &t - edges;
	462	}/}}}/
[12365]	463	/FUNCTION Geometry::GetId(const GeomEdge t){{{*/
[5573]	464	long Geometry::GetId(const GeomEdge * t) const {
[5400]	465	return t - edges;
	466	}/}}}/
[12365]	467	/FUNCTION Geometry::GetId(const Curve c){{{*/
[5400]	468	long Geometry::GetId(const Curve * c) const {
	469	return c - curves;
	470	}/}}}/
[12365]	471	/FUNCTION Geometry::Containing{{{/
[5573]	472	GeomEdge* Geometry::Containing(const R2 P, GeomEdge * start) const {
[5400]	473	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)/
	474
[5573]	475	GeomEdge* on =start,* pon=0;
[5400]	476	// walk with the cos on geometry
	477	int counter=0;
	478	while(pon != on){
	479	counter++;
[6412]	480	_assert_(counter<100);
[5400]	481	pon = on;
	482	R2 A= (*on)[0];
	483	R2 B= (*on)[1];
	484	R2 AB = B-A;
	485	R2 AP = P-A;
	486	R2 BP = P-B;
	487	if ( (AB,AP) < 0)
	488	on = on->Adj[0];
	489	else if ( (AB,BP) > 0)
	490	on = on->Adj[1];
	491	else
	492	return on;
	493	}
	494	return on;
	495	}
[12365]	496	/}}}/
	497	/FUNCTION Geometry::PostRead{{{/
[5397]	498	void Geometry::PostRead(){
[3913]	499	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)/
	500
[5177]	501	long i,j,k;
	502	int jj;
	503	long *head_v = new long[nbv];
	504	long next_p = new long[2nbe];
	505	float *eangle = new float[nbe];
	506	double eps = 1e-20;
[12218]	507	BamgQuadtree quadtree; // build quadtree to find duplicates
[5177]	508	BamgVertex *v0 = vertices;
[5573]	509	GeomVertex v0g = (GeomVertex) (void*)v0;
[3913]	510
	511	k=0;
	512
	513	//build quadtree for this geometry
	514	for (i=0;i<nbv;i++){
	515
	516	/*build integer coordinates (non unique)
[5400]	517	these coordinates are used by the quadtree to group
	518	the vertices by groups of 5:
	519	All the coordinates are transformed to ]0,1[^2
	520	then, the integer coordinates are computed using
	521	the transformation ]0,1[^2 -> [0 2^30-1[^2 for a quadtree of depth 30*/
[5180]	522	vertices[i].i=R2ToI2(vertices[i].r);
[3913]	523
[5177]	524	/find nearest vertex already present in the quadtree (NULL if empty)/
[5120]	525	BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);
[3913]	526
[5177]	527	/if there is a vertex found that is to close to vertices[i] -> error/
	528	if( v && Norme1(v->r - vertices[i].r) < eps ){
[12509]	529	_printLine_("reference numbers: " << v->ReferenceNumber << " " << vertices[i].ReferenceNumber);
	530	_printLine_("Id: " << i+1);
[13385]	531	_printLine_("Coords: ["<<v->r.x<<" "<<v->r.y<<"] ["<<vertices[i].r.x<<" "<<vertices[i].r.y<<"]");
	532
[9306]	533	delete [] next_p;
	534	delete [] head_v;
	535	delete [] eangle;
[13036]	536	_error_("two points of the geometry are very closed to each other (see reference numbers above)");
[3913]	537	}
	538
[5177]	539	/Add vertices[i] to the quadtree/
[5124]	540	quadtree.Add(vertices[i]);
[3913]	541	}
	542
	543	/* Here we use a powerful chaining algorithm
	544	*
	545	* 1. What is a chaining algorithm?
	546	*
	547	* If F is a function that goes from i in [0 n] to j in [0 m]
	548	* and we want to compute the reciprocal function F-1 of F
	549	* (what are the antecedents of a given j in Im(F) )
	550	* We use 2 lists:
	551	* head_F[j] that holds the head of lists
	552	* next_F[i] that holds the list of elements that have the same image
	553	*
	554	* Example:
	555	* i1, i2, ..., ip in [0,n] are all antecedents of a given j in [0 m]
	556	* head_F[j] = ip
	557	* next_F[ip]= ip-1
	558	* ....
	559	* next_F[i2]= i1
	560	* next_F[i1]= -1 //end of the list
	561	*
	562	* Algorithm:
	563	* for(j=0;j<m;j++) head_F[j] = -1 //initialization
	564	* for(i=0;i<n;i++){
	565	* j=F[i];
	566	* next_F[i]= head_F[j];
	567	* head_F[j]=i;
	568	* }
	569	*
	570	* Then, we can go through all the elements that have for image j:
	571	* for(i=head_F[j]; i!=-1; i=next_F[i])
	572	* initialization of i by i=head_F[j]
	573	* stop the loop when i=-1 (end of the chain)
	574	* iterate using i=next_F[i] (next element that have for image j)
	575	*
	576	* 2. How to use this algorithm here?
	577	*
	578	* Here F is a function that associates two vertices v0 and v1 for a given edge E
	579	* We want to build the reciprocal function: what are the edges that contains
	580	* a vertex v?
	581	* To do so, we use the same chaining algorithm but there is a difficulty
[5339]	582	* coming from the fact that for F we have a couple of vertices and not one
	583	* vertex.
	584	* To overcome this difficulty, we use a global indexing exactly like in
[3913]	585	* C/C++ so that
	586	* a member of a 2-column-table can be described by one index p=i*2+j
	587	* i=(int)p/2 line number of p
[5339]	588	* j=p%2 column number of p
[3913]	589	*
	590	* Algorithm:
	591	* for(i=0;i<nbv;i++) head_v[i] = -1 //initialization
	592	* for(i=0;i<nbe;i++){
	593	* for(j=0;j<2;j++){
	594	* p=2*i+j;
	595	* v=edges(i,j);
	596	* next_p[p]= head_v[v];
	597	* head_v[v]=p;
	598	* }
	599	* }
	600	*/
	601
	602	//initialize head_v as -1
	603	for (i=0;i<nbv;i++) head_v[i]=-1;
	604	k=0;
	605	for (i=0;i<nbe;i++) {
	606	//compute vector of edge i that goes from vertex 0 to vertex 1
	607	R2 v10=edges[i].v[1]->r - edges[i].v[0]->r;
	608	double lv10=Norme2(v10);
	609	//check that its length is not 0
[9309]	610	if(lv10==0){
	611	delete [] next_p;
	612	delete [] head_v;
	613	delete [] eangle;
[13036]	614	_error_("Length of edge " << i << " is 0");
[9309]	615	}
[3913]	616	//compute angle in [-Pi Pi]
	617	eangle[i] = atan2(v10.y,v10.x);
	618	//build chains head_v and next_p
	619	for (j=0;j<2;j++){
[5149]	620	long v=GetId(edges[i].v[j]);
[3913]	621	next_p[k]=head_v[v];
	622	head_v[v]=k++; //post increment: head_v[v]=k; and then k=k+1;
	623	}
	624	}
	625
	626	//sort head_v by order of increasing edges angle
	627	for (i=0;i<nbv;i++) {
[5339]	628	int exch=1,ord=0;
[3913]	629
	630	//exchange vertices position in head_v and next_p till tey are sorted
	631	while (exch){
[5339]	632	long *p=head_v+i;
	633	long *po=p;
	634	long n=*p;
[3913]	635	register float angleold=-1000 ; // angle = - infinity
	636	ord=0; exch=0;
	637
[5339]	638	// loop over the edges that contain the vertex i (till -1)
[3913]	639	while (n >=0){
	640	ord++;
[5339]	641	long i1=n/2; // i1 = floor (n/2) -> Edge number
	642	long j1=n%2; // j1 = 1 if n is odd -> Vertex index for this edge (0 or 1)
	643	long* pn=next_p+n;
[3913]	644
[5339]	645	//Next vertex index
[3913]	646	n=*pn;
	647
	648	//compute angle between horizontal axis and v0->v1
	649	float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1];
	650
	651	//exchange if the current edge angle is smaller than the previous one
	652	if (angleold > angle){
	653	exch=1;
	654	pn=po; // next_p[n] = n + 1
	655	po=p; //
	656	*p=n; // next_p[n+1] = n
	657	po=pn; // po now point toward pn (invert next and current)
	658	}
	659
	660	//else, continue going to the next edge position
	661	else{ // to have : po -> p -> pn
	662	angleold=angle; // update maximum angle
	663	po=p; // po now point toward p (current position)
	664	p=pn; // p now point toward pn (next position)
	665	}
	666	}
	667	}
	668
[5339]	669	// angular test on current vertex to guess whether it is a corner (ord = number of edges holding i)
	670	if(ord==2) {
[3913]	671	long n1 = head_v[i];
	672	long n2 = next_p[n1];
	673	long i1 = n1/2, i2 = n2/2; // edge number
	674	long j1 = n1%2, j2 = n2%2; // vertex in the edge
	675	float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1];
	676	float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2];
	677	float da12 = Abs(angle2-angle1);
[5091]	678	if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) {
[3913]	679	vertices[i].SetCorner() ;
	680	}
[5339]	681	// if the edge type/referencenumber a changing then is SetRequired();
[5146]	682	if (edges[i1].type != edges[i2].type \|\| edges[i1].Required()){
[3913]	683	vertices[i].SetRequired();
	684	}
[5148]	685	if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) {
[3913]	686	vertices[i].SetRequired();
	687	}
	688	}
	689	if(ord != 2) {
	690	vertices[i].SetCorner();
	691	}
	692
[5339]	693	/close the list around the vertex to have a circular loop/
[3913]	694	long no=-1, ne = head_v[i];
	695	while (ne >=0) ne = next_p[no=ne];
	696	if(no>=0) next_p[no] = head_v[i];
	697	}
	698
[5339]	699	/*Check that the list makes sense (we have all the time the same vertex)
	700	* and build adjacent edges*/
[3913]	701	k =0;
	702	for (i=0;i<nbe;i++){
	703	for (j=0;j<2;j++){
[5339]	704
[3913]	705	long n1 = next_p[k++];
	706	long i1 = n1/2 ,j1=n1%2;
[5339]	707
[13036]	708	if( edges[i1].v[j1] != edges[i].v[j]) _error_("Problem while processing edges: check the edge list");
[5339]	709
[3913]	710	edges[i1].Adj[j1] = edges + i;
[5401]	711	edges[i1].AdjVertexIndex[j1] = j;
[3913]	712	}
	713	}
	714
[5397]	715	/* generation of all the tangents*/
[3913]	716	for (i=0;i<nbe;i++) {
[5339]	717	R2 AB =edges[i].v[1]->r -edges[i].v[0]->r;// AB = vertex0 -> vertex1
	718	double lAB=Norme2(AB); // Get length of AB
	719	double ltg2[2]={0.0}; // initialize tangent
[3913]	720
	721	//loop over the 2 vertices of the edge
[5339]	722	for (j=0;j<2;j++) {
	723	R2 tg =edges[i].tg[j];
[3913]	724	double ltg=Norme2(tg);
	725
	726	//by default, tangent=[0 0]
[5339]	727	if(ltg==0){
[3913]	728	//if the current vertex of the edge is not a corner
[5339]	729	if(!edges[i].v[j]->Corner()){
	730	/*The tangent is set as the vector between the
	731	* previous and next vertices connected to current vertex
	732	* normed by the edge length*/
[5401]	733	tg = edges[i].v[1-j]->r - edges[i].Adj[j]->v[1-edges[i].AdjVertexIndex[j]]->r;
[3913]	734	ltg= Norme2(tg);
	735	tg = tg *(lAB/ltg);
	736	ltg= lAB;
	737	}
[5339]	738	//else: a Corner no tangent => nothing to do
[3913]	739	}
	740	else{
	741	//tangent has already been computed
	742	tg = tg*(lAB/ltg),ltg=lAB;
	743	}
	744
[5339]	745	ltg2[j] = ltg;
[3913]	746
[5339]	747	if ((tg,AB)<0) tg = -tg;
[3913]	748
[5339]	749	edges[i].tg[j]=tg;
[3913]	750	}
	751	if (ltg2[0]!=0) edges[i].SetTgA();
	752	if (ltg2[1]!=0) edges[i].SetTgB();
	753	}
	754
[5397]	755	/* generation of all curves (from corner to corner)*/
	756	/We proceed in 2 steps: first allocate, second build/
[3913]	757	for (int step=0;step<2;step++){
[5339]	758
[5397]	759	//unmark all edges
[3913]	760	for (i=0;i<nbe;i++) edges[i].SetUnMark();
[5397]	761	long nb_marked_edges=0;
[5339]	762
[5397]	763	//initialize number of curves
[5148]	764	nbcurves = 0;
[4997]	765
[5397]	766	for (int level=0;level<2 && nb_marked_edges!=nbe;level++){
[5339]	767	for (i=0;i<nbe;i++){
[5397]	768
[5573]	769	GeomEdge & ei=edges[i];
[5339]	770	for(j=0;j<2;j++){
[5397]	771	/*If current edge ei is unmarked and (level=1 or vertex i is required (corner)):
	772	* we do have the first edge of a new curve*/
[5339]	773	if (!ei.Mark() && (level \|\| ei[j].Required())) {
	774	int k0=j,k1;
[5573]	775	GeomEdge *e=&ei;
	776	GeomVertex a=(e)(k0); // begin
[5339]	777	if(curves){
[5400]	778	curves[nbcurves].FirstEdge=e;
	779	curves[nbcurves].FirstVertexIndex=k0;
[5339]	780	}
	781	int nee=0;
[5397]	782	for(;;){
[5339]	783	nee++;
	784	k1 = 1-k0; // next vertex of the edge
	785	e->SetMark();
[5397]	786	nb_marked_edges++;
[5339]	787	e->CurveNumber=nbcurves;
[5573]	788	GeomVertex b=(e)(k1);
[5397]	789
	790	//break if we have reached the other end of the curve
	791	if (a==b \|\| b->Required()){
[5400]	792	if(curves){
	793	curves[nbcurves].LastEdge=e;
	794	curves[nbcurves].LastVertexIndex=k1;
	795	}
[5397]	796	break;
	797	}
	798	//else: go to next edge (adjacent)
	799	else{
[5401]	800	k0 = e->AdjVertexIndex[k1];// vertex in next edge
[5397]	801	e = e->Adj[k1]; // next edge
	802	}
[5339]	803	}
	804	nbcurves++;
	805	if(level) a->SetRequired();
	806	}
	807	}
	808	}
	809	}
[6412]	810	_assert_(nb_marked_edges && nbe);
[5397]	811	//allocate if first step
	812	if(step==0) curves=new Curve[nbcurves];
[3913]	813	}
	814
	815	/clean up/
[4997]	816	delete [] next_p;
	817	delete [] head_v;
	818	delete [] eangle;
[3913]	819
	820	}
[12365]	821	/}}}/
	822	/FUNCTION Geometry::ProjectOnCurve {{{/
[5573]	823	GeomEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
[3913]	824	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)/
	825	/Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge/
	826
	827	double save_s=s;
	828	int NbTry=0;
	829
	830	retry:
	831
	832	s=save_s;
[5573]	833	GeomEdge* on=e.GeomEdgeHook;
[3913]	834	if (!on){
[13036]	835	_error_("ProjectOnCurve error message: edge provided should be on geometry");
[3913]	836	}
[5573]	837	if (!e[0].GeomEdgeHook \|\| !e[1].GeomEdgeHook){
[13036]	838	_error_("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry");
[3913]	839	}
	840
	841	//Get the two vertices of the edge
[5120]	842	const BamgVertex &v0=e[0];
	843	const BamgVertex &v1=e[1];
[3913]	844
	845	//Get position of V0, V1 and vector v0->v1
	846	R2 V0=v0,V1=v1,V01=V1-V0;
	847
	848	//Get geometrical vertices corresponding to v0 and v1
[5573]	849	VertexOnGeom vg0=v0.GeomEdgeHook, vg1=v1.GeomEdgeHook;
[3913]	850
	851	//build two pointers towrad current geometrical edge
[5573]	852	GeomEdge eg0=on, eg1=on;
[3913]	853
	854	//Get edge direction and swap v0 and v1 if necessary
	855	R2 Ag=(R2)(on)[0],Bg=(R2)(on)[1],AB=Bg-Ag;
	856	int OppositeSens = (V01,AB)<0;
[5148]	857	int direction0=0,direction1=1;
[3913]	858	if (OppositeSens) s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);
	859
	860	//Compute metric of new vertex as a linear interpolation of the two others
	861	V.m=Metric(1.0-s,v0,s,v1);
	862
	863	const int mxe=100;
[5573]	864	GeomEdge* ge[mxe+1];
[5148]	865	int directionge[mxe+1];
[3913]	866	double lge[mxe+1];
	867	int bge=mxe/2,tge=bge;
[5573]	868	ge[bge] = e.GeomEdgeHook;
[5148]	869	directionge[bge]=1;
[3913]	870
[5573]	871	while (eg0!=(GeomEdge)vg0 && (eg0)(direction0)!=(GeomVertex*)vg0){
[3913]	872	if (bge<=0) {
	873	if(NbTry) {
[12509]	874	_printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
	875	_printLine_("That bug might come from:");
	876	_printLine_(" 1) a mesh edge containing more than " << mxe/2 << " geometrical edges");
	877	_printLine_(" 2) code bug : be sure that we call Mesh::SetVertexFieldOn() before");
	878	_printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
[13036]	879	_error_("see above");
[5130]	880	}
[3913]	881	NbTry++;
	882	goto retry;
	883	}
[5573]	884	GeomEdge* tmpge = eg0;
[5148]	885	ge[--bge] =eg0 = eg0->Adj[direction0];
[6412]	886	_assert_(bge>=0 && bge<=mxe);
[5401]	887	direction0 = 1-( directionge[bge] = tmpge->AdjVertexIndex[direction0]);
[3913]	888	}
[5573]	889	while (eg1 != (GeomEdge) vg1 && (eg1)(direction1) != (GeomVertex*) vg1) {
[3913]	890	if(tge>=mxe ) {
[12509]	891	_printLine_("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)");
[3913]	892	NbTry++;
	893	if (NbTry<2) goto retry;
[12509]	894	_printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
	895	_printLine_("That bug might come from:");
	896	_printLine_(" 1) a mesh edge contening more than " << mxe/2 << " geometrical edges");
	897	_printLine_(" 2) code bug : be sure that we call Mesh::SetVertexFieldOn() before");
	898	_printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
[13036]	899	_error_("see above");
[3913]	900	}
[5573]	901	GeomEdge* tmpge = eg1;
[5148]	902	ge[++tge] =eg1 = eg1->Adj[direction1];
[5401]	903	directionge[tge]= direction1 = 1-tmpge->AdjVertexIndex[direction1];
[6412]	904	_assert_(tge>=0 && tge<=mxe);
[3913]	905	}
	906
	907
[5573]	908	if ((eg0)(direction0)==(GeomVertex)vg0)
[5148]	909	vg0=VertexOnGeom((BamgVertex) vg0,*eg0,direction0); //vg0 = absisce
[3913]	910
[5573]	911	if ((eg1)(direction1)==(GeomVertex)vg1)
[5148]	912	vg1=VertexOnGeom((BamgVertex) vg1,*eg1,direction1);
[3913]	913
	914	double sg;
	915	if (eg0 == eg1) {
	916	register double s0=vg0,s1=vg1;
	917	sg = s0(1.0-s) + ss1;
	918	on=eg0;
	919	}
	920	else {
	921	R2 AA=V0,BB;
	922	double s0,s1;
	923	int i=bge;
	924	double ll=0;
	925	for(i=bge;i<tge;i++){
[6412]	926	_assert_(i>=0 && i<=mxe);
[5148]	927	BB = (*ge[i])[directionge[i]];
[3913]	928	lge[i]=ll += Norme2(AA-BB);
	929	AA=BB ;}
	930	lge[tge]=ll+=Norme2(AA-V1);
	931	// search the geometrical edge
[6412]	932	_assert_(s<=1.0);
[3913]	933	double ls= s*ll;
	934	on =0;
	935	s0 = vg0;
[5148]	936	s1= directionge[bge];
[3913]	937	double l0=0,l1;
	938	i=bge;
	939	while ( (l1=lge[i]) < ls ) {
[6412]	940	_assert_(i>=0 && i<=mxe);
[5148]	941	i++,s0=1-(s1=directionge[i]),l0=l1;
[3913]	942	}
	943	on=ge[i];
	944	if (i==tge)
	945	s1=vg1;
	946
	947	s =(ls-l0)/(l1-l0);
	948	sg =s0(1.0-s)+ss1;
	949	}
[6412]	950	_assert_(on);
[3913]	951	V.r= on->F(sg);
	952	GV=VertexOnGeom(V,*on,sg);
	953	return on;
	954	}
[12365]	955	/}}}/
	956	/FUNCTION Geometry::R2ToI2{{{/
[5180]	957	I2 Geometry::R2ToI2(const R2 & P) const {
[5177]	958	/*coefIcoor is the coefficient used for integer coordinates:
	959	* (x-pmin.x)
	960	* Icoor x = (2^30 -1) ------------
	961	* D
	962	* where D is the longest side of the domain (direction x or y)
	963	* so that (x-pmin.x)/D is in ]0 1[
	964	*
	965	* coefIcoor = (2^30 -1)/D
	966	*/
	967	return I2( (Icoor1) (coefIcoor(P.x-pmin.x)) ,(Icoor1) (coefIcoor(P.y-pmin.y)) );
[5130]	968	}/}}}/
[12365]	969	/FUNCTION Geometry::UnMarkEdges{{{/
[5130]	970	void Geometry::UnMarkEdges() {
	971	for (int i=0;i<nbe;i++) edges[i].SetUnMark();
	972	}/}}}/
[3913]	973	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: