/[escript]/trunk/finley/src/Mesh_findMatchingFaces.c
ViewVC logotype

Diff of /trunk/finley/src/Mesh_findMatchingFaces.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2747 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 2748 by gross, Tue Nov 17 07:32:59 2009 UTC
# Line 22  Line 22 
22    
23  #include "Util.h"  #include "Util.h"
24  #include "Mesh.h"  #include "Mesh.h"
25    
26  /**************************************************************/  /**************************************************************/
27    
28  static double  Finley_Mesh_lockingGridSize=0;  static double  Finley_Mesh_lockingGridSize=0;
# Line 56  void Finley_Mesh_findMatchingFaces(Finle Line 57  void Finley_Mesh_findMatchingFaces(Finle
57        _dist_=0; \        _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)])); \        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];      char error_msg[LenErrorMsg_MAX];
62      double h=DBLE(HUGE_VAL),h_local,dist,*X=NULL;      double h=DBLE(HUGE_VAL),h_local,dist,*X=NULL;  
     dim_t NN=faces->ReferenceElement->Type->numNodes;  
     dim_t numDim=nodes->numDim;  
63      Finley_Mesh_findMatchingFaces_center *center;      Finley_Mesh_findMatchingFaces_center *center;
64      index_t e_0,e_1,*a1=NULL,*a2=NULL,*perm=NULL,*perm_tmp=NULL,*itmp_ptr=NULL;      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;      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(NN*numDim*faces->numElements,double);      X=TMPMEMALLOC(NN*numDim*faces->numElements,double);
83      center=TMPMEMALLOC(faces->numElements,Finley_Mesh_findMatchingFaces_center);      center=TMPMEMALLOC(faces->numElements,Finley_Mesh_findMatchingFaces_center);
84      a1=TMPMEMALLOC(NN,int);      a1=TMPMEMALLOC(NN,int);
# Line 77  void Finley_Mesh_findMatchingFaces(Finle Line 91  void Finley_Mesh_findMatchingFaces(Finle
91              /* get the element center */              /* get the element center */
92              center[e].refId=e;              center[e].refId=e;
93              for (i=0;i<MAX_numDim;i++) center[e].x[i]=0;              for (i=0;i<MAX_numDim;i++) center[e].x[i]=0;
94              for (i0=0;i0<faces->ReferenceElement->Type->numNodesOnFace;i0++) {              for (i0=0;i0<numNodesOnFace;i0++) {
95                 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]+=X[INDEX3(i,faceNodes[i0],e,numDim,NN)];
96              }              }
97              for (i=0;i<numDim;i++) center[e].x[i]/=faces->ReferenceElement->Type->numNodesOnFace;              for (i=0;i<numDim;i++) center[e].x[i]/=numNodesOnFace;
98              /* get the minimum distance between nodes in the element */              /* get the minimum distance between nodes in the element */
99              for (i0=0;i0<faces->ReferenceElement->Type->numNodesOnFace;i0++) {              for (i0=0;i0<numNodesOnFace;i0++) {
100                 for (i1=i0+1;i1<faces->ReferenceElement->Type->numNodesOnFace;i1++) {                 for (i1=i0+1;i1<numNodesOnFace;i1++) {
101                    getDist(h_local,e,faces->ReferenceElement->Type->faceNode[i0],e,faces->ReferenceElement->Type->faceNode[i1]);                    getDist(h_local,e,faceNodes[i0],e,faceNodes[i1]);
102                    h=MIN(h,h_local);                    h=MIN(h,h_local);
103                 }                 }
104              }              }
# Line 116  void Finley_Mesh_findMatchingFaces(Finle Line 130  void Finley_Mesh_findMatchingFaces(Finle
130                   /* if node 0 and perm[0] are the same we are ready */                   /* if node 0 and perm[0] are the same we are ready */
131                   getDist(dist,e_0,0,e_1,perm[0]);                   getDist(dist,e_0,0,e_1,perm[0]);
132                   if (dist<=h*tolerance) break;                   if (dist<=h*tolerance) break;
133                   if (faces->ReferenceElement->Type->shiftNodes[0]>=0) {                   if (shiftNodes[0]>=0) {
134                      /* rotate the nodes */                      /* rotate the nodes */
135                      itmp_ptr=perm;                      itmp_ptr=perm;
136                      perm=perm_tmp;                      perm=perm_tmp;
137                      perm_tmp=itmp_ptr;                      perm_tmp=itmp_ptr;
138                      #pragma ivdep                      #pragma ivdep
139                      for (i=0;i<NN;i++) perm[i]=perm_tmp[faces->ReferenceElement->Type->shiftNodes[i]];                      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: */                   /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
142                   if (perm[0]==0) {                   if (perm[0]==0) {
# Line 132  void Finley_Mesh_findMatchingFaces(Finle Line 146  void Finley_Mesh_findMatchingFaces(Finle
146                }                }
147                /* now we check if the second nodes match */                /* now we check if the second nodes match */
148                if (Finley_noError()) {                if (Finley_noError()) {
149                   if (faces->ReferenceElement->Type->numNodesOnFace>1) {                   if (numNodesOnFace>1) {
150                      getDist(dist,e_0,1,e_1,perm[faces->ReferenceElement->Type->faceNode[1]]);                      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 */                      /* if the second node does not match we reverse the direction of the nodes */
152                      if (dist>h*tolerance) {                      if (dist>h*tolerance) {
153                            /* rotate the nodes */                            /* rotate the nodes */
154                            if (faces->ReferenceElement->Type->reverseNodes[0]<0) {                            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);                               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);                               Finley_setError(VALUE_ERROR,error_msg);
157                            } else {                            } else {
# Line 145  void Finley_Mesh_findMatchingFaces(Finle Line 159  void Finley_Mesh_findMatchingFaces(Finle
159                               perm=perm_tmp;                               perm=perm_tmp;
160                               perm_tmp=itmp_ptr;                               perm_tmp=itmp_ptr;
161                               #pragma ivdep                               #pragma ivdep
162                               for (i=0;i<NN;i++) perm[i]=perm_tmp[faces->ReferenceElement->Type->reverseNodes[i]];                               for (i=0;i<NN;i++) perm[i]=perm_tmp[reverseNodes[i]];
163                               getDist(dist,e_0,1,e_1,perm[faces->ReferenceElement->Type->faceNode[1]]);                               getDist(dist,e_0,1,e_1,perm[faceNodes[1]]);
164                               if (dist>h*tolerance) {                               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);                                   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);                                   Finley_setError(VALUE_ERROR,error_msg);
# Line 157  void Finley_Mesh_findMatchingFaces(Finle Line 171  void Finley_Mesh_findMatchingFaces(Finle
171                }                }
172                /* we check if the rest of the face nodes match: */                /* we check if the rest of the face nodes match: */
173                if (Finley_noError()) {                if (Finley_noError()) {
174                   for (i=2;i<faces->ReferenceElement->Type->numNodesOnFace;i++) {                   for (i=2;i<numNodesOnFace;i++) {
175                      n=faces->ReferenceElement->Type->faceNode[i];                      n=faceNodes[i];
176                      getDist(dist,e_0,n,e_1,perm[n]);                      getDist(dist,e_0,n,e_1,perm[n]);
177                      if (dist>h*tolerance) {                      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);                         sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the %d-th node of element %d to touching element %d",i,e_0,e_1);

Legend:
Removed from v.2747  
changed lines
  Added in v.2748

  ViewVC Help
Powered by ViewVC 1.1.26