finley/src/Mesh_findMatchingFaces.c


/*******************************************************
*
* Copyright (c) 2003-2009 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


/**************************************************************/

/*   Finley: Mesh */

/* searches for faces in the mesh which are matching */

/**************************************************************/

#include "Util.h"
#include "Mesh.h"

/**************************************************************/

static double  Finley_Mesh_lockingGridSize=0;

int Finley_Mesh_findMatchingFaces_compar(const void *arg1 , const void *arg2 ) {
   Finley_Mesh_findMatchingFaces_center *e1,*e2;
   bool_t l,g;
   dim_t i;
   e1=(Finley_Mesh_findMatchingFaces_center*) arg1;
   e2=(Finley_Mesh_findMatchingFaces_center*) arg2;
   for (i=0;i<MAX_numDim;i++) {
     l= (e1->x[i]<e2->x[i]+Finley_Mesh_lockingGridSize) ? TRUE : FALSE;
     g= (e2->x[i]<e1->x[i]+Finley_Mesh_lockingGridSize) ? TRUE : FALSE;
     if (! (l && g)) {
        if (l) return -1;
        if (g) return 1;
     }
   }
   if ( e1->refId < e2->refId ) {
       return -1;
   } else if ( e1->refId > e2->refId ) {
       return 1;
   } else {
       return 0;
   }
}

void Finley_Mesh_findMatchingFaces(Finley_NodeFile *nodes, Finley_ElementFile *faces, double safety_factor,double tolerance,
                                   dim_t* numPairs, index_t* elem0,index_t* elem1,index_t* matching_nodes_in_elem1) {
#define getDist(_dist_,_e0_,_i0_,_e1_,_i1_) \
      {dim_t i;   \
      _dist_=0; \
      for (i=0;i<numDim;i++) _dist_=MAX(_dist_,ABS(X[INDEX3(i,_i0_,_e0_,numDim,NN)]-X[INDEX3(i,_i1_,_e1_,numDim,NN)])); \
      } 
    Finley_ReferenceElement*  refElement=NULL;
    char error_msg[LenErrorMsg_MAX];
    double h=DBLE(HUGE_VAL),h_local,dist,*X=NULL;   
    Finley_Mesh_findMatchingFaces_center *center;
    index_t e_0,e_1,*a1=NULL,*a2=NULL,*perm=NULL,*perm_tmp=NULL,*itmp_ptr=NULL, *faceNodes=NULL, *shiftNodes=NULL, *reverseNodes=NULL ;
    dim_t e,i,i0,i1,n,NN,numNodesOnFace;
    
    dim_t numDim=nodes->numDim;

    refElement= Finley_ReferenceElementSet_borrowReferenceElement(faces->referenceElementSet, FALSE);
    NN=faces->numNodes;
    
    numNodesOnFace=refElement->Type->numNodesOnFace;
    faceNodes=refElement->Type->faceNodes;
    shiftNodes=refElement->Type->shiftNodes;
    reverseNodes=refElement->Type->reverseNodes;
    
    if (numNodesOnFace<=0) {
     sprintf(error_msg,"Finley_Mesh_findMatchingFaces: matching faces cannot be applied to face elements of type %s",refElement->Type->Name);
     Finley_setError(TYPE_ERROR,error_msg);
     return;
    }
    X=TMPMEMALLOC(NN*numDim*faces->numElements,double);
    center=TMPMEMALLOC(faces->numElements,Finley_Mesh_findMatchingFaces_center);
    a1=TMPMEMALLOC(NN,int);
    a2=TMPMEMALLOC(NN,int);
    if (!(Finley_checkPtr(X) || Finley_checkPtr(center) || Finley_checkPtr(a1) || Finley_checkPtr(a2)) ) {
       /* OMP */
       for (e=0;e<faces->numElements;e++) {
            /* get the coordinates of the nodes */
            Finley_Util_Gather_double(NN,&(faces->Nodes[INDEX2(0,e,NN)]),numDim,nodes->Coordinates,&(X[INDEX3(0,0,e,numDim,NN)]));
            /* get the element center */
            center[e].refId=e;
            for (i=0;i<MAX_numDim;i++) center[e].x[i]=0;
            for (i0=0;i0<numNodesOnFace;i0++) {
               for (i=0;i<numDim;i++) center[e].x[i]+=X[INDEX3(i,faceNodes[i0],e,numDim,NN)];
            }
            for (i=0;i<numDim;i++) center[e].x[i]/=numNodesOnFace;
            /* get the minimum distance between nodes in the element */
            for (i0=0;i0<numNodesOnFace;i0++) {
               for (i1=i0+1;i1<numNodesOnFace;i1++) {
                  getDist(h_local,e,faceNodes[i0],e,faceNodes[i1]);
                  h=MIN(h,h_local);
               }
            }
       }
       /* set the */
       Finley_Mesh_lockingGridSize=h*MAX(safety_factor,0);
       #ifdef Finley_TRACE
       printf("locking grid size is %e\n",Finley_Mesh_lockingGridSize);
       printf("absolute tolerance is %e.\n",h * tolerance);
       #endif
       /* sort the elements by center center coordinates (lexigraphical)*/
       qsort(center,faces->numElements,sizeof(Finley_Mesh_findMatchingFaces_center),Finley_Mesh_findMatchingFaces_compar);
       /* find elements with matching center */
       *numPairs=0;
       /* OMP */
       for (e=0;e<faces->numElements-1 && Finley_noError();e++) {
          dist=0; 
          for (i=0;i<numDim;i++) dist=MAX(dist,ABS(center[e].x[i]-center[e+1].x[i]));
          if (dist < h * tolerance) { 
              e_0=center[e].refId;
              e_1=center[e+1].refId;
              elem0[*numPairs]=e_0;
              elem1[*numPairs]=e_1;
              /* now the element e_1 is rotated such that the first node in element e_0 and e_1 have the same coordinates */
              perm=a1;
              perm_tmp=a2;
              for (i=0;i<NN;i++) perm[i]=i;
              while (Finley_noError()) {
                 /* if node 0 and perm[0] are the same we are ready */
                 getDist(dist,e_0,0,e_1,perm[0]);
                 if (dist<=h*tolerance) break;
                 if (shiftNodes[0]>=0) {
                    /* rotate the nodes */
                    itmp_ptr=perm;
                    perm=perm_tmp;
                    perm_tmp=itmp_ptr;
                    #pragma ivdep
                    for (i=0;i<NN;i++) perm[i]=perm_tmp[shiftNodes[i]];
                 }
                 /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
                 if (perm[0]==0) {
                       sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match first node of element %d to touching element %d",e_0,e_1);
                       Finley_setError(VALUE_ERROR,error_msg);
                 }
              }
              /* now we check if the second nodes match */
              if (Finley_noError()) {
                 if (numNodesOnFace>1) {
                    getDist(dist,e_0,1,e_1,perm[faceNodes[1]]);
                    /* if the second node does not match we reverse the direction of the nodes */
                    if (dist>h*tolerance) {
                          /* rotate the nodes */
                          if (reverseNodes[0]<0) {
                             sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the second node of element %d to touching element %d",e_0,e_1);
                             Finley_setError(VALUE_ERROR,error_msg);
                          } else {
                             itmp_ptr=perm;
                             perm=perm_tmp;
                             perm_tmp=itmp_ptr;
                             #pragma ivdep
                             for (i=0;i<NN;i++) perm[i]=perm_tmp[reverseNodes[i]];
                             getDist(dist,e_0,1,e_1,perm[faceNodes[1]]);
                             if (dist>h*tolerance) {
                                 sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the second node of element %d to touching element %d",e_0,e_1);
                                 Finley_setError(VALUE_ERROR,error_msg);
                             }
                          }
                    }
                 }
              }
              /* we check if the rest of the face nodes match: */
              if (Finley_noError()) {
                 for (i=2;i<numNodesOnFace;i++) {
                    n=faceNodes[i];
                    getDist(dist,e_0,n,e_1,perm[n]);
                    if (dist>h*tolerance) {
                       sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the %d-th node of element %d to touching element %d",i,e_0,e_1);
                       Finley_setError(VALUE_ERROR,error_msg);
                       break;
                    }
                 }
              }
              /* copy over the permuted nodes of e_1 into matching_nodes_in_elem1 */
              if (Finley_noError()) {
                 for (i=0;i<NN;i++)  matching_nodes_in_elem1[INDEX2(i,*numPairs,NN)]=faces->Nodes[INDEX2(perm[i],e_1,NN)];
              }
              (*numPairs)++;
          }
       }
       #ifdef Finley_TRACE
       printf("number of pairs of matching faces %d\n",*numPairs);
       #endif
    }
    /* clean up */
    TMPMEMFREE(X);
    TMPMEMFREE(center);
    TMPMEMFREE(a1);
    TMPMEMFREE(a1);

#undef getDist
}
1
2	/*******************************************************
3	*
4	* Copyright (c) 2003-2009 by University of Queensland
5	* Earth Systems Science Computational Center (ESSCC)
6	* http://www.uq.edu.au/esscc
7	*
8	* Primary Business: Queensland, Australia
9	* Licensed under the Open Software License version 3.0
10	* http://www.opensource.org/licenses/osl-3.0.php
11	*
12	*******************************************************/
13
14
15	/**************************************************************/
16
17	/* Finley: Mesh */
18
19	/* searches for faces in the mesh which are matching */
20
21	/**************************************************************/
22
23	#include "Util.h"
24	#include "Mesh.h"
25
26	/**************************************************************/
27
28	static double Finley_Mesh_lockingGridSize=0;
29
30	int Finley_Mesh_findMatchingFaces_compar(const void arg1 , const void arg2 ) {
31	Finley_Mesh_findMatchingFaces_center e1,e2;
32	bool_t l,g;
33	dim_t i;
34	e1=(Finley_Mesh_findMatchingFaces_center*) arg1;
35	e2=(Finley_Mesh_findMatchingFaces_center*) arg2;
36	for (i=0;i<MAX_numDim;i++) {
37	l= (e1->x[i]<e2->x[i]+Finley_Mesh_lockingGridSize) ? TRUE : FALSE;
38	g= (e2->x[i]<e1->x[i]+Finley_Mesh_lockingGridSize) ? TRUE : FALSE;
39	if (! (l && g)) {
40	if (l) return -1;
41	if (g) return 1;
42	}
43	}
44	if ( e1->refId < e2->refId ) {
45	return -1;
46	} else if ( e1->refId > e2->refId ) {
47	return 1;
48	} else {
49	return 0;
50	}
51	}
52
53	void Finley_Mesh_findMatchingFaces(Finley_NodeFile nodes, Finley_ElementFile faces, double safety_factor,double tolerance,
54	dim_t* numPairs, index_t* elem0,index_t* elem1,index_t* matching_nodes_in_elem1) {
55	#define getDist(_dist_,_e0_,_i0_,_e1_,_i1_) \
56	{dim_t i; \
57	_dist_=0; \
58	for (i=0;i<numDim;i++) _dist_=MAX(_dist_,ABS(X[INDEX3(i,_i0_,_e0_,numDim,NN)]-X[INDEX3(i,_i1_,_e1_,numDim,NN)])); \
59	}
60	Finley_ReferenceElement* refElement=NULL;
61	char error_msg[LenErrorMsg_MAX];
62	double h=DBLE(HUGE_VAL),h_local,dist,*X=NULL;
63	Finley_Mesh_findMatchingFaces_center *center;
64	index_t e_0,e_1,a1=NULL,a2=NULL,perm=NULL,perm_tmp=NULL,itmp_ptr=NULL, faceNodes=NULL, shiftNodes=NULL, reverseNodes=NULL ;
65	dim_t e,i,i0,i1,n,NN,numNodesOnFace;
66
67	dim_t numDim=nodes->numDim;
68
69	refElement= Finley_ReferenceElementSet_borrowReferenceElement(faces->referenceElementSet, FALSE);
70	NN=faces->numNodes;
71
72	numNodesOnFace=refElement->Type->numNodesOnFace;
73	faceNodes=refElement->Type->faceNodes;
74	shiftNodes=refElement->Type->shiftNodes;
75	reverseNodes=refElement->Type->reverseNodes;
76
77	if (numNodesOnFace<=0) {
78	sprintf(error_msg,"Finley_Mesh_findMatchingFaces: matching faces cannot be applied to face elements of type %s",refElement->Type->Name);
79	Finley_setError(TYPE_ERROR,error_msg);
80	return;
81	}
82	X=TMPMEMALLOC(NNnumDimfaces->numElements,double);
83	center=TMPMEMALLOC(faces->numElements,Finley_Mesh_findMatchingFaces_center);
84	a1=TMPMEMALLOC(NN,int);
85	a2=TMPMEMALLOC(NN,int);
86	if (!(Finley_checkPtr(X) \|\| Finley_checkPtr(center) \|\| Finley_checkPtr(a1) \|\| Finley_checkPtr(a2)) ) {
87	/* OMP */
88	for (e=0;e<faces->numElements;e++) {
89	/* get the coordinates of the nodes */
90	Finley_Util_Gather_double(NN,&(faces->Nodes[INDEX2(0,e,NN)]),numDim,nodes->Coordinates,&(X[INDEX3(0,0,e,numDim,NN)]));
91	/* get the element center */
92	center[e].refId=e;
93	for (i=0;i<MAX_numDim;i++) center[e].x[i]=0;
94	for (i0=0;i0<numNodesOnFace;i0++) {
95	for (i=0;i<numDim;i++) center[e].x[i]+=X[INDEX3(i,faceNodes[i0],e,numDim,NN)];
96	}
97	for (i=0;i<numDim;i++) center[e].x[i]/=numNodesOnFace;
98	/* get the minimum distance between nodes in the element */
99	for (i0=0;i0<numNodesOnFace;i0++) {
100	for (i1=i0+1;i1<numNodesOnFace;i1++) {
101	getDist(h_local,e,faceNodes[i0],e,faceNodes[i1]);
102	h=MIN(h,h_local);
103	}
104	}
105	}
106	/* set the */
107	Finley_Mesh_lockingGridSize=h*MAX(safety_factor,0);
108	#ifdef Finley_TRACE
109	printf("locking grid size is %e\n",Finley_Mesh_lockingGridSize);
110	printf("absolute tolerance is %e.\n",h * tolerance);
111	#endif
112	/* sort the elements by center center coordinates (lexigraphical)*/
113	qsort(center,faces->numElements,sizeof(Finley_Mesh_findMatchingFaces_center),Finley_Mesh_findMatchingFaces_compar);
114	/* find elements with matching center */
115	*numPairs=0;
116	/* OMP */
117	for (e=0;e<faces->numElements-1 && Finley_noError();e++) {
118	dist=0;
119	for (i=0;i<numDim;i++) dist=MAX(dist,ABS(center[e].x[i]-center[e+1].x[i]));
120	if (dist < h * tolerance) {
121	e_0=center[e].refId;
122	e_1=center[e+1].refId;
123	elem0[*numPairs]=e_0;
124	elem1[*numPairs]=e_1;
125	/* now the element e_1 is rotated such that the first node in element e_0 and e_1 have the same coordinates */
126	perm=a1;
127	perm_tmp=a2;
128	for (i=0;i<NN;i++) perm[i]=i;
129	while (Finley_noError()) {
130	/* if node 0 and perm[0] are the same we are ready */
131	getDist(dist,e_0,0,e_1,perm[0]);
132	if (dist<=h*tolerance) break;
133	if (shiftNodes[0]>=0) {
134	/* rotate the nodes */
135	itmp_ptr=perm;
136	perm=perm_tmp;
137	perm_tmp=itmp_ptr;
138	#pragma ivdep
139	for (i=0;i<NN;i++) perm[i]=perm_tmp[shiftNodes[i]];
140	}
141	/* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
142	if (perm[0]==0) {
143	sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match first node of element %d to touching element %d",e_0,e_1);
144	Finley_setError(VALUE_ERROR,error_msg);
145	}
146	}
147	/* now we check if the second nodes match */
148	if (Finley_noError()) {
149	if (numNodesOnFace>1) {
150	getDist(dist,e_0,1,e_1,perm[faceNodes[1]]);
151	/* if the second node does not match we reverse the direction of the nodes */
152	if (dist>h*tolerance) {
153	/* rotate the nodes */
154	if (reverseNodes[0]<0) {
155	sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the second node of element %d to touching element %d",e_0,e_1);
156	Finley_setError(VALUE_ERROR,error_msg);
157	} else {
158	itmp_ptr=perm;
159	perm=perm_tmp;
160	perm_tmp=itmp_ptr;
161	#pragma ivdep
162	for (i=0;i<NN;i++) perm[i]=perm_tmp[reverseNodes[i]];
163	getDist(dist,e_0,1,e_1,perm[faceNodes[1]]);
164	if (dist>h*tolerance) {
165	sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the second node of element %d to touching element %d",e_0,e_1);
166	Finley_setError(VALUE_ERROR,error_msg);
167	}
168	}
169	}
170	}
171	}
172	/* we check if the rest of the face nodes match: */
173	if (Finley_noError()) {
174	for (i=2;i<numNodesOnFace;i++) {
175	n=faceNodes[i];
176	getDist(dist,e_0,n,e_1,perm[n]);
177	if (dist>h*tolerance) {
178	sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the %d-th node of element %d to touching element %d",i,e_0,e_1);
179	Finley_setError(VALUE_ERROR,error_msg);
180	break;
181	}
182	}
183	}
184	/* copy over the permuted nodes of e_1 into matching_nodes_in_elem1 */
185	if (Finley_noError()) {
186	for (i=0;i<NN;i++) matching_nodes_in_elem1[INDEX2(i,*numPairs,NN)]=faces->Nodes[INDEX2(perm[i],e_1,NN)];
187	}
188	(*numPairs)++;
189	}
190	}
191	#ifdef Finley_TRACE
192	printf("number of pairs of matching faces %d\n",*numPairs);
193	#endif
194	}
195	/* clean up */
196	TMPMEMFREE(X);
197	TMPMEMFREE(center);
198	TMPMEMFREE(a1);
199	TMPMEMFREE(a1);
200
201	#undef getDist
202	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision