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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1062 - (show annotations)
Mon Mar 26 06:17:53 2007 UTC (12 years, 6 months ago) by gross
File MIME type: text/plain
File size: 6660 byte(s)
reduced integration schemes are implemented now for grad, integrate, etc. Tests still to be added.
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 /* detects faces in the mesh that match and replaces it by step elements */
18
19 /**************************************************************/
20
21 /* Author: gross@access.edu.au */
22 /* Version: $Id$ */
23
24 /**************************************************************/
25
26 #include "Mesh.h"
27
28 /**************************************************************/
29
30 void Finley_Mesh_joinFaces(Finley_Mesh* self,double safety_factor,double tolerance, bool_t optimize_labeling) {
31
32 char error_msg[LenErrorMsg_MAX];
33 index_t e0,e1,*elem1=NULL,*elem0=NULL,*elem_mask=NULL,*matching_nodes_in_elem1=NULL;
34 Finley_ElementFile *newFaceElementsFile=NULL,*newContactElementsFile=NULL;
35 dim_t e,i,numPairs, NN, NN_Contact,c, new_numFaceElements;
36 if (self->FaceElements==NULL) return;
37
38 if (self->FaceElements->ReferenceElement->Type->numNodesOnFace<=0) {
39 sprintf(error_msg,"Finley_Mesh_joinFaces:joining faces cannot be applied to face elements of type %s",self->FaceElements->ReferenceElement->Type->Name);
40 Finley_setError(TYPE_ERROR,error_msg);
41 return;
42 }
43 if (self->ContactElements==NULL) {
44 Finley_setError(TYPE_ERROR,"Finley_Mesh_joinFaces: no contact element file present.");
45 return;
46 }
47
48 NN=self->FaceElements->ReferenceElement->Type->numNodes;
49 NN_Contact=self->ContactElements->ReferenceElement->Type->numNodes;
50
51 if (2*NN!=NN_Contact) {
52 sprintf(error_msg,"Finley_Mesh_joinFaces:contact element file for %s cannot hold elements created from face elements %s",
53 self->ContactElements->ReferenceElement->Type->Name,self->FaceElements->ReferenceElement->Type->Name);
54 Finley_setError(TYPE_ERROR,error_msg);
55 return;
56 }
57
58 /* allocate work arrays */
59 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
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 if (Finley_noError()) {
68 /* 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 new_numFaceElements=0;
76 /* 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 #ifndef PASO_MPI
85 newContactElementsFile=Finley_ElementFile_alloc(self->ContactElements->ReferenceElement->Type->TypeId,self->ContactElements->order,self->ContactElements->reduced_order);
86 newFaceElementsFile=Finley_ElementFile_alloc(self->FaceElements->ReferenceElement->Type->TypeId,self->FaceElements->order,self->FaceElements->reduced_order);
87 #else
88 /* TODO */
89 PASO_MPI_TODO;
90 #endif
91 if (Finley_noError()) {
92 Finley_ElementFile_allocTable(newContactElementsFile,numPairs+self->ContactElements->numElements);
93 Finley_ElementFile_allocTable(newFaceElementsFile,new_numFaceElements);
94 }
95 /* copy the old elements over */
96 if (Finley_noError()) {
97 /* 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 c=self->ContactElements->numElements;
102 /* 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 newContactElementsFile->minColor=0;
114 newContactElementsFile->maxColor=numPairs-1;
115 newContactElementsFile->isPrepared=self->FaceElements->isPrepared;
116 }
117 /* set new face and Contact elements */
118 if (Finley_noError()) {
119
120 Finley_ElementFile_dealloc(self->FaceElements);
121 self->FaceElements=newFaceElementsFile;
122 Finley_ElementFile_prepare(&(self->FaceElements),self->Nodes->numNodes,self->Nodes->degreeOfFreedom);
123
124 Finley_ElementFile_dealloc(self->ContactElements);
125 self->ContactElements=newContactElementsFile;
126 Finley_ElementFile_prepare(&(self->ContactElements),self->Nodes->numNodes,self->Nodes->degreeOfFreedom);
127
128 Finley_Mesh_prepare(self);
129
130 } else {
131 Finley_ElementFile_dealloc(newFaceElementsFile);
132 Finley_ElementFile_dealloc(newContactElementsFile);
133 }
134 }
135 }
136 TMPMEMFREE(elem1);
137 TMPMEMFREE(elem0);
138 TMPMEMFREE(matching_nodes_in_elem1);
139 TMPMEMFREE(elem_mask);
140 if (Finley_noError()) {
141 if ( ! Finley_Mesh_isPrepared(self) ) {
142 Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
143 }
144 }
145 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26