finley/src/Mesh_findMatchingFaces.c

/*
 ************************************************************
 *          Copyright 2006 by ACcESS MNRF                   *
 *                                                          *
 *              http://www.access.edu.au                    *
 *       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 */

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

/*  Author: gross@access.edu.au */
/*  Version: $Id$

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

#include "Util.h"
#include "Mesh.h"
#ifdef USE_QSORTG
#include "qsortG.h"
#endif
/**************************************************************/

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;
     }
   }
   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)])); \
      } 

#define SWAP(_i1_,_i2_) \
            {index_t* i;  \
              i=(_i2_); \
              (_i2_)=(_i1_); \
              (_i1_)=i; \
             }

    char error_msg[LenErrorMsg_MAX];
    double h=DBLE(HUGE_VAL),h_local,dist,*X=NULL;
    dim_t NN=faces->ReferenceElement->Type->numNodes;
    dim_t numDim=nodes->numDim;
    Finley_Mesh_findMatchingFaces_center *center;
    index_t e_0,e_1,a1[NN],a2[NN],*perm,*perm_tmp;
    dim_t e,i,i0,i1,n;

    X=TMPMEMALLOC(NN*numDim*faces->numElements,double);
    center=TMPMEMALLOC(faces->numElements,Finley_Mesh_findMatchingFaces_center);
    if (!(Finley_checkPtr(X) || Finley_checkPtr(center)) ) {
       /* 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<faces->ReferenceElement->Type->numNodesOnFace;i0++) {
               for (i=0;i<numDim;i++) center[e].x[i]+=X[INDEX3(i,faces->ReferenceElement->Type->faceNode[i0],e,numDim,NN)];
            }
            for (i=0;i<numDim;i++) center[e].x[i]/=faces->ReferenceElement->Type->numNodesOnFace;
            /* get the minimum distance between nodes in the element */
            for (i0=0;i0<faces->ReferenceElement->Type->numNodesOnFace;i0++) {
               for (i1=i0+1;i1<faces->ReferenceElement->Type->numNodesOnFace;i1++) {
                  getDist(h_local,e,faces->ReferenceElement->Type->faceNode[i0],e,faces->ReferenceElement->Type->faceNode[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);
       #endif
       /* sort the elements by center center coordinates (lexigraphical)*/
#ifdef USE_QSORTG
       qsortG(center,faces->numElements,sizeof(Finley_Mesh_findMatchingFaces_center),Finley_Mesh_findMatchingFaces_compar);
#else
       qsort(center,faces->numElements,sizeof(Finley_Mesh_findMatchingFaces_center),Finley_Mesh_findMatchingFaces_compar);
#endif
       /* find elements with matching center */
       *numPairs=0;
       /* OMP */
       for (e=0;e<faces->numElements-1 && Finley_noError();e++) {
          if (Finley_Mesh_findMatchingFaces_compar((void*) &(center[e]),(void*) &(center[e+1]))==0) {
              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 (faces->ReferenceElement->Type->shiftNodes[0]>=0) {
                    /* rotate the nodes */
                    for (i=0;i<NN;i++) perm_tmp[i]=perm[faces->ReferenceElement->Type->shiftNodes[i]];
                    SWAP(perm,perm_tmp);
                 }
                 /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
                 if (perm[0]==0) {
                       sprintf(error_msg,"__FILE__: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 (faces->ReferenceElement->Type->numNodesOnFace>1) {
                    getDist(dist,e_0,1,e_1,perm[faces->ReferenceElement->Type->faceNode[1]]);
                    /* if the second node does not match we reverse the direction of the nodes */
                    if (dist>h*tolerance) {
                          /* rotate the nodes */
                          if (faces->ReferenceElement->Type->reverseNodes[0]<0) {
                             sprintf(error_msg,"__FILE__:couldn't match the second node of element %d to touching element %d",e_0,e_1);
                             Finley_setError(VALUE_ERROR,error_msg);
                          } else {
                             for (i=0;i<NN;i++) perm_tmp[i]=perm[faces->ReferenceElement->Type->reverseNodes[i]];
                             SWAP(perm,perm_tmp);
                             getDist(dist,e_0,1,e_1,perm[faces->ReferenceElement->Type->faceNode[1]]);
                             if (dist>h*tolerance) {
                                 sprintf(error_msg,"__FILE__: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<faces->ReferenceElement->Type->numNodesOnFace;i++) {
                    n=faces->ReferenceElement->Type->faceNode[i];
                    getDist(dist,e_0,n,e_1,perm[n]);
                    if (dist>h*tolerance) {
                       sprintf(error_msg,"__FILE__: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);

#undef getDist
#undef SWAP
}

/*
* $Log$
* Revision 1.6  2005/09/15 03:44:22  jgs
* Merge of development branch dev-02 back to main trunk on 2005-09-15
*
* Revision 1.5.2.1  2005/09/07 06:26:19  gross
* the solver from finley are put into the standalone package paso now
*
* Revision 1.5  2005/07/08 04:07:51  jgs
* Merge of development branch back to main trunk on 2005-07-08
*
* Revision 1.4  2004/12/15 07:08:32  jgs
* *** empty log message ***
* Revision 1.1.1.1.2.2  2005/06/29 02:34:51  gross
* some changes towards 64 integers in finley
*
* Revision 1.1.1.1.2.1  2004/11/24 01:37:14  gross
* some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
*
*
*
*/

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