Context Navigation

Mesh2.h@ 3232

Last change on this file since 3232 was 3232, checked in by Mathieu Morlighem, 15 years ago
Separated all objects (still some work to do)
File size: 14.8 KB

Rev	Line
[2865]	1	#ifndef _MESH2_H_
	2	#define _MESH2_H_
	3
	4	#include "../objects/objects.h"
	5	#include "../shared/shared.h"
	6	#include "../include/macros.h"
	7	#include "../toolkits/toolkits.h"
	8
	9	#include "meshtype.h"
	10
[3232]	11	#include "objects/Metric.h"
	12	#include "objects/DoubleAndInt4.h"
	13	#include "objects/Direction.h"
	14	#include "objects/Vertex.h"
	15	#include "objects/TriangleAdjacent.h"
	16	#include "objects/Edge.h"
	17	#include "objects/GeometricalVertex.h"
	18	#include "objects/GeometricalEdge.h"
	19	#include "objects/Curve.h"
	20	#include "objects/Triangle.h"
	21	#include "objects/ListofIntersectionTriangles.h"
	22	#include "objects/GeometricalSubDomain.h"
	23	#include "objects/SubDomain.h"
	24	#include "objects/VertexOnGeom.h"
	25	#include "objects/VertexOnVertex.h"
	26	#include "objects/VertexOnEdge.h"
	27	#include "objects/CrackedEdge.h"
	28	#include "objects/Triangles.h"
	29	#include "objects/Geometry.h"
	30
[2865]	31	namespace bamg {
[2992]	32
	33	/INLINE functions{{{1/
[2870]	34	// to sort in descending order
	35	template<class T> inline void HeapSort(T *c,long n){
[3022]	36	/Intermediary/
	37	int l,j,r,i;
	38	T crit;
	39	c--; //the array must starts at 1 and not 0
	40	if(n<=1) return; //return if size <=1
	41	l=n/2+1; //initialize l and r
	42	r=n;
	43	for(;;){
	44	if(l<=1){
	45	crit =c[r];
	46	c[r--]=c[1];
	47	if (r==1){c[1]=crit; return;}
	48	}
	49	else crit = c[--l];
[2870]	50	j=l;
[3022]	51	for(;;){
[2870]	52	i=j;
	53	j=2*j;
[3022]	54	if (j>r) {c[i]=crit;break;}
	55	if ((j<r) && (c[j] < c[j+1])) j++;//c[j+1]> c[j] -> take j+1 instead of j (larger value)
	56	if (crit < c[j]) c[i]=c[j]; //c[j] > crit -> stock this large value in i(<j)
	57	else{c[i]=crit;break;} //c[j] < crit -> stock crit in i (<j)
[2870]	58	}
	59	}
	60	}
[3022]	61	// to sort in descending order and return ordering
	62	template<class T> inline void HeapSort(int** porder,T* c,int n){
	63	/Intermediary/
	64	int l,j,r,i;
	65	T crit;
	66	int pos;
	67	int* order=NULL;
	68	order=(int)xmalloc(nsizeof(int));
	69	for(i=0;i<n;i++) order[i]=i+1;
	70	c--; //the array must starts at 1 and not 0
	71	order--;
	72	if(n<=1) return; //return if size <=1
	73	l=n/2+1; //initialize l and r
	74	r=n;
	75	for(;;){
	76	if(l<=1){
	77	crit =c[r]; pos=order[r];
	78	c[r--]=c[1]; order[r+1]=order[1];
	79	if (r==1){
	80	c[1]=crit; order[1]=pos;
	81	order++;
	82	*porder=order;
	83	return;
	84	}
	85	}
	86	else {crit=c[--l]; pos=order[l];}
	87	j=l;
	88	for(;;){
	89	i=j;
	90	j=2*j;
	91	if (j>r) {c[i]=crit;order[i]=pos;break;}
	92	if ((j<r) && (c[j] < c[j+1]))j++;
	93	if (crit < c[j]) {c[i]=c[j];order[i]=order[j];}
	94	else{c[i]=crit;order[i]=pos;break;}
	95	}
	96	}
	97	}
[2983]	98	inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3]){
	99	return A[0]( B[1]C[2]-B[2]*C[1])
[2992]	100	- A[1]( B[0]C[2]-B[2]*C[0])
	101	+ A[2]( B[0]C[1]-B[1]*C[0]);
[2983]	102	}
[2992]	103	inline TriangleAdjacent Adj(const TriangleAdjacent & a)
	104	{ return a.Adj();}
[2865]	105
[2992]	106	inline TriangleAdjacent Next(const TriangleAdjacent & ta)
	107	{ return TriangleAdjacent(ta.t,NextEdge[ta.a]);}
	108
	109	inline TriangleAdjacent Previous(const TriangleAdjacent & ta)
	110	{ return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);}
	111	inline void Adj(GeometricalEdge * & on,int &i)
	112	{int j=i;i=on->DirAdj[i];on=on->Adj[j];}
	113	inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc)
	114	{
	115	Real4 ret;
	116	I2 ia=va,ib=vb,ic=vc;
	117	I2 ab=ib-ia,bc=ic-ib,ac=ic-ia;
	118	Icoor2 deta=Det(ab,ac);
	119	if (deta <=0) ret = -1;
	120	else {
	121	Real8 a = sqrt((Real8) (ac,ac)),
	122	b = sqrt((Real8) (bc,bc)),
	123	c = sqrt((Real8) (ab,ab)),
	124	p = a+b+c;
	125	Real8 h= Max(Max(a,b),c),ro=deta/p;
	126	ret = ro/h;}
	127	return ret;
	128	}
	129	Icoor2 inline det(const Vertex & a,const Vertex & b,const Vertex & c){
	130	register Icoor2 bax = b.i.x - a.i.x ,bay = b.i.y - a.i.y;
	131	register Icoor2 cax = c.i.x - a.i.x ,cay = c.i.y - a.i.y;
	132	return baxcay - baycax;
	133	}
	134	inline TriangleAdjacent FindTriangleAdjacent(Edge &E){
	135	Vertex * a = E.v[0];
	136	Vertex * b = E.v[1];
	137
	138	Triangle * t = a->t;
	139	int i = a->vint;
	140	TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge
	141	if (!t \|\| i<0 \|\| i>=3){
	142	throw ErrorException(__FUNCT__,exprintf("!t \|\| i<0 !! i>=3"));
	143	}
	144	if ( a!=(*t)(i)){
	145	throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)"));
	146	}
	147	int k=0;
	148	do { // turn around vertex in direct sens (trigo)
	149	k++;
	150	if (k>=20000){
	151	throw ErrorException(__FUNCT__,exprintf("k>=20000"));
	152	}
	153	// in no crack => ta.EdgeVertex(1) == a otherwise ???
	154	if (ta.EdgeVertex(1) == a && ta.EdgeVertex(0) == b) return ta; // find
	155	ta = ta.Adj();
	156	if (ta.EdgeVertex(0) == a && ta.EdgeVertex(1) == b) return ta; // find
	157	--ta;
	158	} while (t != (Triangle *)ta);
	159	throw ErrorException(__FUNCT__,exprintf("FindTriangleAdjacent: triangle not found"));
	160	return TriangleAdjacent(0,0);//for compiler
	161	}
	162	inline Vertex* TheVertex(Vertex * a){
	163	// give a unique vertex with smallest number
	164	// in case on crack in mesh
	165	Vertex * r(a), *rr;
	166	Triangle * t = a->t;
	167	int i = a->vint;
	168	TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge
	169	if (!t \|\| i<0 \|\| i>=3){
	170	throw ErrorException(__FUNCT__,exprintf("!t \|\| i<0 !! i>=3"));
	171	}
	172	if ( a!=(*t)(i)){
	173	throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)"));
	174	}
	175	int k=0;
	176	do { // turn around vertex in direct sens (trigo)
	177	k++;
	178	if (k>=20000){
	179	throw ErrorException(__FUNCT__,exprintf("k>=20000"));
	180	}
	181	// in no crack => ta.EdgeVertex(1) == a
	182	if ((rr=ta.EdgeVertex(0)) < r) r = rr;
	183	ta = ta.Adj();
	184	if ((rr=ta.EdgeVertex(1)) < r) r =rr;
	185	--ta;
	186	} while (t != (Triangle*) ta);
	187	return r;
	188	}
	189
	190	/}}}1/
	191
	192	/INLINE functions of CLASS GeometricalVertex{{{1/
	193	inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ){
	194	*this = rec;
	195	}
	196	/}}}1/
	197	/INLINE functions of CLASS VertexOnVertex{{{1/
	198	inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) {
	199	*this = Th.VertexOnBThVertex[i];
	200	v = ThNew.vertices + Th.Number(v);
	201
	202	}
	203	/}}}1/
	204	/INLINE functions of CLASS Triangles{{{1/
[3022]	205	inline void Triangles::ReMakeTriangleContainingTheVertex(){
[2992]	206	register Int4 i;
[3022]	207	for ( i=0;i<nbv;i++){
	208	vertices[i].vint=0;
[2992]	209	vertices[i].t=0;
	210	}
[3022]	211	for ( i=0;i<nbt;i++) triangles[i].SetTriangleContainingTheVertex();
[2992]	212	}
[2865]	213
[2992]	214	inline void Triangles::UnMarkUnSwapTriangle()
	215	{
	216	register Int4 i;
	217	for ( i=0;i<nbt;i++)
	218	for(int j=0;j<3;j++)
	219	triangles[i].SetUnMarkUnSwap(j);
	220	}
[2865]	221
[2992]	222	inline void Triangles::SetVertexFieldOn(){
	223	for (Int4 i=0;i<nbv;i++)
	224	vertices[i].onGeometry=0;
	225	for (Int4 j=0;j<NbVerticesOnGeomVertex;j++ )
	226	VerticesOnGeomVertex[j].SetOn();
	227	for (Int4 k=0;k<NbVerticesOnGeomEdge;k++ )
	228	VerticesOnGeomEdge[k].SetOn();
	229	}
	230	inline void Triangles::SetVertexFieldOnBTh(){
	231	for (Int4 i=0;i<nbv;i++)
	232	vertices[i].onGeometry=0;
	233	for (Int4 j=0;j<NbVertexOnBThVertex;j++ )
	234	VertexOnBThVertex[j].SetOnBTh();
	235	for (Int4 k=0;k<NbVertexOnBThEdge;k++ )
	236	VertexOnBThEdge[k].SetOnBTh();
[2865]	237
[2992]	238	}
	239	/}}}1/
	240	/INLINE functions of CLASS Triangle{{{1/
[2870]	241	inline Triangle* Triangle::Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const
	242	{
	243	// return the other triangle of the quad if a quad or 0 if not a quat
	244	Triangle * t =0;
	245	if (link) {
	246	int a=-1;
[2970]	247	if (TriaAdjSharedEdge[0] & 16 ) a=0;
	248	if (TriaAdjSharedEdge[1] & 16 ) a=1;
	249	if (TriaAdjSharedEdge[2] & 16 ) a=2;
[2870]	250	if (a>=0) {
[2970]	251	t = TriaAdjTriangles[a];
[2870]	252	// if (t-this<0) return 0;
[2970]	253	v2 = TriaVertices[VerticesOfTriangularEdge[a][0]];
	254	v0 = TriaVertices[VerticesOfTriangularEdge[a][1]];
	255	v1 = TriaVertices[OppositeEdge[a]];
	256	v3 = t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
[2870]	257	}
	258	}
	259	return t;
	260	}
[2865]	261
[2870]	262	inline double Triangle::QualityQuad(int a,int option) const
	263	{ // first do the logique part
	264	double q;
[2970]	265	if (!link \|\| TriaAdjSharedEdge[a] &4)
[2870]	266	q= -1;
	267	else {
[2970]	268	Triangle * t = TriaAdjTriangles[a];
[2870]	269	if (t-this<0) q= -1;// because we do 2 times
	270	else if (!t->link ) q= -1;
[2970]	271	else if (TriaAdjSharedEdge[0] & 16 \|\| TriaAdjSharedEdge[1] & 16 \|\| TriaAdjSharedEdge[2] & 16 \|\| t->TriaAdjSharedEdge[0] & 16 \|\| t->TriaAdjSharedEdge[1] & 16 \|\| t->TriaAdjSharedEdge[2] & 16 )
[2870]	272	q= -1;
	273	else if(option)
	274	{
[2970]	275	const Vertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];
	276	const Vertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];
	277	const Vertex & v1 = *TriaVertices[OppositeEdge[a]];
	278	const Vertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
[2870]	279	q = QuadQuality(v0,v1,v2,v3); // do the float part
	280	}
	281	else q= 1;
	282	}
	283	return q;
	284	}
[2992]	285	inline void Triangle::Set(const Triangle & rec,const Triangles & Th ,Triangles & ThNew)
	286	{
	287	*this = rec;
	288	if ( TriaVertices[0] ) TriaVertices[0] = ThNew.vertices + Th.Number(TriaVertices[0]);
	289	if ( TriaVertices[1] ) TriaVertices[1] = ThNew.vertices + Th.Number(TriaVertices[1]);
	290	if ( TriaVertices[2] ) TriaVertices[2] = ThNew.vertices + Th.Number(TriaVertices[2]);
	291	if(TriaAdjTriangles[0]) TriaAdjTriangles[0] = ThNew.triangles + Th.Number(TriaAdjTriangles[0]);
	292	if(TriaAdjTriangles[1]) TriaAdjTriangles[1] = ThNew.triangles + Th.Number(TriaAdjTriangles[1]);
	293	if(TriaAdjTriangles[2]) TriaAdjTriangles[2] = ThNew.triangles + Th.Number(TriaAdjTriangles[2]);
	294	if (link >= Th.triangles && link < Th.triangles + Th.nbt)
	295	link = ThNew.triangles + Th.Number(link);
	296	}
	297	inline Triangle::Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k) {
	298	Vertex *v=Th->vertices;
	299	Int4 nbv = Th->nbv;
	300	if (i<0 \|\| j<0 \|\| k<0){
	301	throw ErrorException(__FUNCT__,exprintf("i<0 \|\| j<0 \|\| k<0"));
	302	}
	303	if (i>=nbv \|\| j>=nbv \|\| k>=nbv){
	304	throw ErrorException(__FUNCT__,exprintf("i>=nbv \|\| j>=nbv \|\| k>=nbv"));
	305	}
	306	TriaVertices[0]=v+i;
	307	TriaVertices[1]=v+j;
	308	TriaVertices[2]=v+k;
	309	TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
	310	TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
	311	det=0;
	312	}
	313	inline Triangle::Triangle(Vertex v0,Vertex v1,Vertex *v2){
	314	TriaVertices[0]=v0;
	315	TriaVertices[1]=v1;
	316	TriaVertices[2]=v2;
	317	TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
	318	TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
	319	if (v0) det=0;
	320	else {
	321	det=-1;
	322	link=NULL;};
	323	}
	324	inline Real4 Triangle::qualite()
	325	{
	326	return det < 0 ? -1 : bamg::qualite(TriaVertices[0],TriaVertices[1],*TriaVertices[2]);
	327	}
[2865]	328
[2992]	329	/}}}1/
	330	/INLINE functions of CLASS Edge{{{1/
[2870]	331	inline void Edge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
	332	{
	333	*this = Th.edges[i];
	334	v[0] = ThNew.vertices + Th.Number(v[0]);
	335	v[1] = ThNew.vertices + Th.Number(v[1]);
[2945]	336	if (onGeometry)
	337	onGeometry = ThNew.Gh.edges+Th.Gh.Number(onGeometry);
[2870]	338	if (adj[0]) adj[0] = ThNew.edges + Th.Number(adj[0]);
	339	if (adj[1]) adj[1] = ThNew.edges + Th.Number(adj[1]);
[2865]	340
[2870]	341	}
[2992]	342
	343	/}}}1/
	344	/INLINE functions of CLASS GeometricalEdge{{{1/
[2870]	345	inline void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew)
	346	{
	347	*this = rec;
	348	v[0] = GhNew.vertices + Gh.Number(v[0]);
	349	v[1] = GhNew.vertices + Gh.Number(v[1]);
	350	if (Adj[0]) Adj[0] = GhNew.edges + Gh.Number(Adj[0]);
	351	if (Adj[1]) Adj[1] = GhNew.edges + Gh.Number(Adj[1]);
	352	}
[2992]	353	/}}}1/
	354	/INLINE functions of CLASS Curve{{{1/
[2870]	355	inline void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew)
	356	{
	357	*this = rec;
	358	be = GhNew.edges + Gh.Number(be);
	359	ee = GhNew.edges + Gh.Number(ee);
	360	if(next) next= GhNew.curves + Gh.Number(next);
	361	}
[2992]	362	/}}}1/
	363	/INLINE functions of CLASS SubDomain{{{1/
[2870]	364	inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
	365	{
	366	*this = Th.subdomains[i];
	367	if ( head-Th.triangles<0 \|\| head-Th.triangles>=Th.nbt){
	368	throw ErrorException(__FUNCT__,exprintf("head-Th.triangles<0 \|\| head-Th.triangles>=Th.nbt"));
	369	}
	370	head = ThNew.triangles + Th.Number(head) ;
	371	if (edge-Th.edges<0 \|\| edge-Th.edges>=Th.nbe);{
	372	throw ErrorException(__FUNCT__,exprintf("edge-Th.edges<0 \|\| edge-Th.edges>=Th.nbe"));
	373	}
	374	edge = ThNew.edges+ Th.Number(edge);
	375	}
[2992]	376	/}}}1/
	377	/INLINE functions of CLASS GeometricalSubDomain{{{1/
[2870]	378	inline void GeometricalSubDomain::Set(const GeometricalSubDomain & rec,const Geometry & Gh ,const Geometry & GhNew)
	379	{
	380	*this = rec;
	381	edge = Gh.Number(edge) + GhNew.edges;
	382	}
[2992]	383	/}}}1/
[3232]	384	/INLINE functions of CLASS Vertex{{{1/
	385	Int4 inline Vertex::Optim(int i,int koption){
	386	Int4 ret=0;
	387	if ( t && (vint >= 0 ) && (vint <3) ){
	388	ret = t->Optim(vint,koption);
	389	if(!i){
	390	t =0; // for no future optime
	391	vint= 0;
	392	}
	393	}
	394	return ret;
	395	}
	396	inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ){
	397	*this = rec;
	398	}
	399	/}}}1/
[2992]	400	/INLINE functions of CLASS VertexOnEdge{{{1/
[2870]	401	inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
	402	{
	403	*this = Th.VertexOnBThEdge[i];
	404	v = ThNew.vertices + Th.Number(v);
	405	}
[2992]	406	/}}}1/
	407	/INLINE functions of CLASS VertexOnGeom{{{1/
[2870]	408	inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew)
	409	{
	410	*this = rec;
	411	mv = ThNew.vertices + Th.Number(mv);
	412	if (gv)
	413	if (abscisse < 0 )
	414	gv = ThNew.Gh.vertices + Th.Gh.Number(gv);
	415	else
	416	ge = ThNew.Gh.edges + Th.Gh.Number(ge);
[2865]	417
[2870]	418	}
[2992]	419	/}}}1/
	420	/INLINE functions of CLASS TriangleAdjacent{{{1/
[2958]	421	inline void TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l ){
	422	//set Adjacent Triangle of a triangle
[2870]	423	if(t) {
[2970]	424	t->TriaAdjTriangles[a]=ta.t;
	425	t->TriaAdjSharedEdge[a]=ta.a\|l;
[2958]	426	}
	427	if(ta.t) {
[2970]	428	ta.t->TriaAdjTriangles[ta.a] = t ;
	429	ta.t->TriaAdjSharedEdge[ta.a] = a\| l ;
[2958]	430	}
	431	}
[2870]	432	inline int TriangleAdjacent::Locked() const
[2970]	433	{ return t->TriaAdjSharedEdge[a] &4;}
[2870]	434	inline int TriangleAdjacent::Cracked() const
[2970]	435	{ return t->TriaAdjSharedEdge[a] &32;}
[2870]	436	inline int TriangleAdjacent::GetAllFlag_UnSwap() const
[2970]	437	{ return t->TriaAdjSharedEdge[a] & 1012;} // take all flag except MarkUnSwap
[2870]	438	inline int TriangleAdjacent::MarkUnSwap() const
[2970]	439	{ return t->TriaAdjSharedEdge[a] &8;}
[2870]	440	inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);}
	441	inline void TriangleAdjacent::SetCracked() { t->SetCracked(a);}
	442	inline TriangleAdjacent TriangleAdjacent::Adj() const
	443	{ return t->Adj(a);}
	444	inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const
[2970]	445	{return t->TriaVertices[VerticesOfTriangularEdge[a][i]]; }
[2870]	446	inline Vertex * TriangleAdjacent::OppositeVertex() const
[2970]	447	{return t->TriaVertices[bamg::OppositeVertex[a]]; }
[2870]	448	inline Icoor2 & TriangleAdjacent::det() const
	449	{ return t->det;}
[2992]	450	int inline TriangleAdjacent::swap()
	451	{ return t->swap(a);}
	452	/}}}1/
[2870]	453
[2992]	454	/Other prototypes {{{1/
	455	TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ;
	456	TriangleAdjacent CloseBoundaryEdgeV2(I2 A,Triangle *t, double &a,double &b);
	457	Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside);
	458	void swap(Triangle *t1,Int1 a1,
	459	Triangle *t2,Int1 a2,
	460	Vertex s1,Vertex s2,Icoor2 det1,Icoor2 det2);
	461	int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
	462	TriangleAdjacent & tt1,Icoor2 & dets1,
	463	Icoor2 & detsa,Icoor2 & detsb, int & nbswap);
	464	int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ;
	465	/}}}1/
[2870]	466
[2865]	467	}
	468	#endif

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

Context Navigation

source: issm/trunk/src/c/Bamgx/Mesh2.h@ 3232

Download in other formats: