/[escript]/branches/doubleplusgood/finley/src/Mesh_joinFaces.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/finley/src/Mesh_joinFaces.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4327 - (show annotations)
Wed Mar 20 05:09:11 2013 UTC (6 years, 6 months ago) by jfenwick
File size: 6193 byte(s)
some finley memory
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Finley: Mesh */
20
21 /* detects faces in the mesh that match and replaces it by step elements */
22
23 /************************************************************************************/
24
25 #include "Mesh.h"
26
27 /************************************************************************************/
28
29 void Finley_Mesh_joinFaces(Finley_Mesh* self,double safety_factor,double tolerance, bool_t optimize) {
30
31 char error_msg[LenErrorMsg_MAX];
32 index_t e0,e1,*elem1=NULL,*elem0=NULL,*elem_mask=NULL,*matching_nodes_in_elem1=NULL;
33 Finley_ElementFile *newFaceElementsFile=NULL,*newContactElementsFile=NULL;
34 dim_t e,i,numPairs, NN, NN_Contact,c, new_numFaceElements;
35 Finley_ReferenceElement* faceRefElement=NULL, *contactRefElement=NULL;
36
37 if (self->MPIInfo->size>1) {
38 Finley_setError(TYPE_ERROR,"Finley_Mesh_joinFaces: MPI is not supported yet.");
39 return;
40 }
41 if (self->ContactElements==NULL) {
42 Finley_setError(TYPE_ERROR,"Finley_Mesh_joinFaces: no contact element file present.");
43 return;
44 }
45 if (self->FaceElements==NULL) return;
46 faceRefElement= Finley_ReferenceElementSet_borrowReferenceElement(self->FaceElements->referenceElementSet, FALSE);
47 contactRefElement= Finley_ReferenceElementSet_borrowReferenceElement(self->ContactElements->referenceElementSet, FALSE);
48
49
50 NN=self->FaceElements->numNodes;
51 NN_Contact=self->ContactElements->numNodes;
52
53 if (faceRefElement->Type->numNodesOnFace<=0) {
54 sprintf(error_msg,"Finley_Mesh_joinFaces: joining faces cannot be applied to face elements of type %s",faceRefElement->Type->Name);
55 Finley_setError(TYPE_ERROR,error_msg);
56 return;
57 }
58
59
60 if (contactRefElement->Type->numNodes != 2*faceRefElement->Type->numNodes) {
61 sprintf(error_msg,"Finley_Mesh_joinFaces: contact element file for %s needs to hold elements created from face elements %s", contactRefElement->Type->Name,faceRefElement->Type->Name);
62 Finley_setError(TYPE_ERROR,error_msg);
63 return;
64 }
65
66 /* allocate work arrays */
67 elem1=new index_t[self->FaceElements->numElements];
68 elem0=new index_t[self->FaceElements->numElements];
69 elem_mask=new index_t[self->FaceElements->numElements];
70 matching_nodes_in_elem1=new index_t[self->FaceElements->numElements*NN];
71
72 if (!(Finley_checkPtr(elem1) || Finley_checkPtr(elem0) || Finley_checkPtr(elem_mask) || Finley_checkPtr(matching_nodes_in_elem1))) {
73
74 /* find the matching face elements */
75 Finley_Mesh_findMatchingFaces(self->Nodes,self->FaceElements,safety_factor,tolerance,&numPairs,elem0,elem1,matching_nodes_in_elem1);
76 if (Finley_noError()) {
77 /* get a list of the face elements to be kept */
78 #pragma omp parallel for private(e) schedule(static)
79 for(e=0;e<self->FaceElements->numElements;e++) elem_mask[e]=1;
80 for(e=0;e<numPairs;e++) {
81 elem_mask[elem0[e]]=0;
82 elem_mask[elem1[e]]=0;
83 }
84 new_numFaceElements=0;
85 /* OMP */
86 for(e=0;e<self->FaceElements->numElements;e++) {
87 if (elem_mask[e]>0) {
88 elem_mask[new_numFaceElements]=e;
89 new_numFaceElements++;
90 }
91 }
92 /* allocate new face element and Contact element files */
93 newContactElementsFile=Finley_ElementFile_alloc(self->ContactElements->referenceElementSet, self->MPIInfo);
94 newFaceElementsFile=Finley_ElementFile_alloc(self->FaceElements->referenceElementSet, self->MPIInfo);
95 if (Finley_noError()) {
96 Finley_ElementFile_allocTable(newContactElementsFile,numPairs+self->ContactElements->numElements);
97 Finley_ElementFile_allocTable(newFaceElementsFile,new_numFaceElements);
98 }
99 /* copy the old elements over */
100 if (Finley_noError()) {
101 /* get the face elements which are still in use:*/
102 Finley_ElementFile_gather(elem_mask,self->FaceElements,newFaceElementsFile);
103 /* get the Contact elements which are still in use:*/
104 Finley_ElementFile_copyTable(0,newContactElementsFile,0,0,self->ContactElements);
105 c=self->ContactElements->numElements;
106 /* OMP */
107 for (e=0;e<numPairs;e++) {
108 e0=elem0[e];
109 e1=elem1[e];
110 newContactElementsFile->Id[c]=MIN(self->FaceElements->Id[e0],self->FaceElements->Id[e1]);
111 newContactElementsFile->Tag[c]=MIN(self->FaceElements->Tag[e0],self->FaceElements->Tag[e1]);
112 newContactElementsFile->Color[c]=e;
113 for (i=0;i<NN;i++) newContactElementsFile->Nodes[INDEX2(i,c,NN_Contact)]=self->FaceElements->Nodes[INDEX2(i,e0,NN)];
114 for (i=0;i<NN;i++) newContactElementsFile->Nodes[INDEX2(i+NN,c,NN_Contact)]=matching_nodes_in_elem1[INDEX2(i,e,NN)];
115 c++;
116 }
117 newContactElementsFile->minColor=0;
118 newContactElementsFile->maxColor=numPairs-1;
119 }
120 /* set new face and Contact elements */
121 if (Finley_noError()) {
122
123 Finley_ElementFile_free(self->FaceElements);
124 self->FaceElements=newFaceElementsFile;
125 Finley_ElementFile_free(self->ContactElements);
126 self->ContactElements=newContactElementsFile;
127 Finley_Mesh_prepare(self, optimize);
128
129 } else {
130 Finley_ElementFile_free(newFaceElementsFile);
131 Finley_ElementFile_free(newContactElementsFile);
132 }
133 }
134 }
135 delete[] elem1;
136 delete[] elem0;
137 delete[] matching_nodes_in_elem1;
138 delete[] elem_mask;
139 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26