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

Annotation of /trunk/finley/src/Mesh_joinFaces.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 730 - (hide annotations)
Mon May 15 04:03:49 2006 UTC (13 years, 6 months ago) by bcumming
File MIME type: text/plain
File size: 6272 byte(s)


1 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12    
13 jgs 82 /**************************************************************/
14    
15     /* Finley: Mesh */
16    
17     /* detects faces in the mesh that match and replaces it by step elements */
18    
19     /**************************************************************/
20    
21 jgs 150 /* Author: gross@access.edu.au */
22     /* Version: $Id$ */
23 jgs 82
24     /**************************************************************/
25    
26     #include "Mesh.h"
27    
28     /**************************************************************/
29    
30     void Finley_Mesh_joinFaces(Finley_Mesh* self,double safety_factor,double tolerance) {
31    
32 jgs 150 char error_msg[LenErrorMsg_MAX];
33 jgs 123 index_t e0,e1,*elem1=NULL,*elem0=NULL,*elem_mask=NULL,*matching_nodes_in_elem1=NULL;
34 jgs 82 Finley_ElementFile *newFaceElementsFile=NULL,*newContactElementsFile=NULL;
35 jgs 123 dim_t e,i,numPairs;
36 jgs 82 if (self->FaceElements==NULL) return;
37    
38     if (self->FaceElements->ReferenceElement->Type->numNodesOnFace<=0) {
39 jgs 153 sprintf(error_msg,"Finley_Mesh_joinFaces:joining faces cannot be applied to face elements of type %s",self->FaceElements->ReferenceElement->Type->Name);
40 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
41 jgs 82 return;
42     }
43     if (self->ContactElements==NULL) {
44 jgs 153 Finley_setError(TYPE_ERROR,"Finley_Mesh_joinFaces: no contact element file present.");
45 jgs 82 return;
46     }
47    
48     int NN=self->FaceElements->ReferenceElement->Type->numNodes;
49     int NN_Contact=self->ContactElements->ReferenceElement->Type->numNodes;
50    
51     if (2*NN!=NN_Contact) {
52 jgs 153 sprintf(error_msg,"Finley_Mesh_joinFaces:contact element file for %s cannot hold elements created from face elements %s",
53 jgs 82 self->ContactElements->ReferenceElement->Type->Name,self->FaceElements->ReferenceElement->Type->Name);
54 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
55 jgs 82 return;
56     }
57    
58     /* allocate work arrays */
59 jgs 123 elem1=TMPMEMALLOC(self->FaceElements->numElements,index_t);
60     elem0=TMPMEMALLOC(self->FaceElements->numElements,index_t);
61     elem_mask=TMPMEMALLOC(self->FaceElements->numElements,index_t);
62     matching_nodes_in_elem1=TMPMEMALLOC(self->FaceElements->numElements*NN,index_t);
63 jgs 82
64     if (!(Finley_checkPtr(elem1) || Finley_checkPtr(elem0) || Finley_checkPtr(elem_mask) || Finley_checkPtr(matching_nodes_in_elem1))) {
65     /* find the matching face elements */
66     Finley_Mesh_findMatchingFaces(self->Nodes,self->FaceElements,safety_factor,tolerance,&numPairs,elem0,elem1,matching_nodes_in_elem1);
67 jgs 150 if (Finley_noError()) {
68 jgs 82 /* get a list of the face elements to be kept */
69     #pragma omp parallel for private(e) schedule(static)
70     for(e=0;e<self->FaceElements->numElements;e++) elem_mask[e]=1;
71     for(e=0;e<numPairs;e++) {
72     elem_mask[elem0[e]]=0;
73     elem_mask[elem1[e]]=0;
74     }
75 jgs 123 dim_t new_numFaceElements=0;
76 jgs 82 /* OMP */
77     for(e=0;e<self->FaceElements->numElements;e++) {
78     if (elem_mask[e]>0) {
79     elem_mask[new_numFaceElements]=e;
80     new_numFaceElements++;
81     }
82     }
83     /* allocate new face element and Contact element files */
84 bcumming 730 #ifndef PASO_MPI
85 jgs 82 newContactElementsFile=Finley_ElementFile_alloc(self->ContactElements->ReferenceElement->Type->TypeId,self->ContactElements->order);
86     newFaceElementsFile=Finley_ElementFile_alloc(self->FaceElements->ReferenceElement->Type->TypeId,self->FaceElements->order);
87 bcumming 730 #else
88     /* TODO */
89     PASO_MPI_TODO;
90     #endif
91 jgs 150 if (Finley_noError()) {
92 jgs 82 Finley_ElementFile_allocTable(newContactElementsFile,numPairs+self->ContactElements->numElements);
93     Finley_ElementFile_allocTable(newFaceElementsFile,new_numFaceElements);
94     }
95     /* copy the old elements over */
96 jgs 150 if (Finley_noError()) {
97 jgs 82 /* get the face elements which are still in use:*/
98     Finley_ElementFile_gather(elem_mask,self->FaceElements,newFaceElementsFile);
99     /* get the Contact elements which are still in use:*/
100     Finley_ElementFile_copyTable(0,newContactElementsFile,0,0,self->ContactElements);
101 jgs 123 dim_t c=self->ContactElements->numElements;
102 jgs 82 /* OMP */
103     for (e=0;e<numPairs;e++) {
104     e0=elem0[e];
105     e1=elem1[e];
106     newContactElementsFile->Id[c]=MIN(self->FaceElements->Id[e0],self->FaceElements->Id[e1]);
107     newContactElementsFile->Tag[c]=MIN(self->FaceElements->Tag[e0],self->FaceElements->Tag[e1]);
108     newContactElementsFile->Color[c]=e;
109     for (i=0;i<NN;i++) newContactElementsFile->Nodes[INDEX2(i,c,NN_Contact)]=self->FaceElements->Nodes[INDEX2(i,e0,NN)];
110     for (i=0;i<NN;i++) newContactElementsFile->Nodes[INDEX2(i+NN,c,NN_Contact)]=matching_nodes_in_elem1[INDEX2(i,e,NN)];
111     c++;
112     }
113 jgs 123 newContactElementsFile->minColor=0;
114     newContactElementsFile->maxColor=numPairs-1;
115 jgs 82 }
116     /* set new face and Contact elements */
117 jgs 150 if (Finley_noError()) {
118 jgs 82
119     Finley_ElementFile_dealloc(self->FaceElements);
120     self->FaceElements=newFaceElementsFile;
121     Finley_ElementFile_prepare(&(self->FaceElements),self->Nodes->numNodes,self->Nodes->degreeOfFreedom);
122    
123     Finley_ElementFile_dealloc(self->ContactElements);
124     self->ContactElements=newContactElementsFile;
125     Finley_ElementFile_prepare(&(self->ContactElements),self->Nodes->numNodes,self->Nodes->degreeOfFreedom);
126    
127 jgs 153 Finley_Mesh_prepare(self);
128 jgs 82
129     } else {
130     Finley_ElementFile_dealloc(newFaceElementsFile);
131     Finley_ElementFile_dealloc(newContactElementsFile);
132     }
133     }
134     }
135     TMPMEMFREE(elem1);
136     TMPMEMFREE(elem0);
137     TMPMEMFREE(matching_nodes_in_elem1);
138     TMPMEMFREE(elem_mask);
139     }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26