/[escript]/branches/schroedinger_upto1946/finley/src/Mesh_joinFaces.c
ViewVC logotype

Contents of /branches/schroedinger_upto1946/finley/src/Mesh_joinFaces.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1947 - (show annotations)
Wed Oct 29 23:19:45 2008 UTC (10 years, 4 months ago) by jfenwick
File MIME type: text/plain
File size: 5968 byte(s)
This does not actually have the changes in it yet.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26