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

Contents of /branches/doubleplusgood/finley/src/Mesh_addPoints.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: 8720 byte(s)
some finley memory
1 /*****************************************************************************
2 *
3 * Copyright (c) 2003-2013 by University of Queensland
4 * http://www.uq.edu.au
5 *
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 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 * Development since 2012 by School of Earth Sciences
12 *
13 *****************************************************************************/
14
15 /************************************************************************************
16
17 Finley: adds points to Points table of the mesh
18
19 ************************************************************************************/
20
21 #include "Mesh.h"
22
23 #ifdef ESYS_MPI
24 void Finley_Mesh_MPI_minimizeDistance( void *, void*, int *, MPI_Datatype * );
25
26 void Finley_Mesh_MPI_minimizeDistance(void *invec_p, void *inoutvec_p, int *len, MPI_Datatype *dtype)
27 {
28 const dim_t numPoints = (*len)/2;
29 double *invec = (double*) invec_p;
30 double *inoutvec = (double*) inoutvec_p;
31 int i;
32 for ( i=0; i<numPoints; i++ ) {
33 if ( invec[2*i] < inoutvec[2*i] ) {
34 inoutvec[2*i]=invec[2*i];
35 inoutvec[2*i+1]=invec[2*i+1];
36 }
37 }
38 }
39 #endif
40
41 void Finley_Mesh_addPoints(Finley_Mesh* mesh, const dim_t numPoints, const double *points_ptr, const index_t *tags_ptr)
42 {
43 dim_t i,n,k, numOldPoints, numNewPoints;
44 double *dist_p=NULL;
45 index_t *node_id_p=NULL, *point_index_p=NULL;
46 index_t order=mesh->integrationOrder;
47 index_t reduced_order=mesh->reducedIntegrationOrder;
48 const dim_t numDim=mesh->Nodes->numDim;
49 Esys_MPIInfo * mpi_info= Esys_MPIInfo_getReference(mesh->MPIInfo);
50 Finley_ElementFile *oldPoints=mesh->Points, *newPoints=NULL;
51 Finley_ReferenceElementSet *refPoints=NULL;
52 const index_t firstDOF=Paso_Distribution_getFirstComponent(mesh->Nodes->degreesOfFreedomDistribution);
53 const index_t lastDOF=Paso_Distribution_getLastComponent(mesh->Nodes->degreesOfFreedomDistribution);
54
55
56 if ( oldPoints == NULL) {
57 refPoints=Finley_ReferenceElementSet_alloc(Finley_Point1, order, reduced_order);
58 numOldPoints=0;
59 } else {
60 refPoints=Finley_ReferenceElementSet_reference(oldPoints->referenceElementSet);
61 numOldPoints=mesh->Points->numElements;
62 }
63 newPoints=Finley_ElementFile_alloc(refPoints, mpi_info);
64 /* first we find the node which is the closest on this processor: */
65 dist_p=new double[numPoints];
66 point_index_p=new index_t[numPoints];
67 node_id_p=new index_t[numPoints];
68
69 for (i=0; i< numPoints; ++i) {
70 dist_p[i]=LARGE_POSITIVE_FLOAT;
71 node_id_p[i]=-1;
72 node_id_p[i]=-1;
73 }
74 if ( numDim == 3 ) {
75 #pragma omp parallel private(i)
76 {
77 for (i=0; i< numPoints; ++i) {
78 register const double X0=points_ptr[INDEX2(0,i,numDim)];
79 register const double X1=points_ptr[INDEX2(1,i,numDim)];
80 register const double X2=points_ptr[INDEX2(2,i,numDim)];
81 register double dist_local=LARGE_POSITIVE_FLOAT;
82 register index_t node_id_local=-1;
83 #pragma omp for private(n)
84 for (n=0; n< mesh-> Nodes->numNodes; n++) {
85 register const double D0=mesh-> Nodes->Coordinates[INDEX2(0,n,numDim)] - X0;
86 register const double D1=mesh-> Nodes->Coordinates[INDEX2(1,n,numDim)] - X1;
87 register const double D2=mesh-> Nodes->Coordinates[INDEX2(2,n,numDim)] - X2;
88 register const double d= D0 * D0 + D1 * D1 + D2 * D2;
89 if ( d < dist_local) {
90 dist_local = d;
91 node_id_local=n;
92 }
93 }
94 #pragma omp critical
95 {
96 if ( dist_local < dist_p[i] ) {
97 dist_p[i] = dist_local;
98 node_id_p[i] = node_id_local;
99 }
100 }
101 }
102 }
103 } else if ( numDim == 2 ) {
104 #pragma omp parallel private(i)
105 {
106 for (i=0; i< numPoints; ++i) {
107 register const double X0=points_ptr[INDEX2(0,i,numDim)];
108 register const double X1=points_ptr[INDEX2(1,i,numDim)];
109 register double dist_local=LARGE_POSITIVE_FLOAT;
110 register index_t node_id_local=-1;
111 #pragma omp for private(n)
112 for (n=0; n< mesh-> Nodes->numNodes; n++) {
113 register const double D0=mesh-> Nodes->Coordinates[INDEX2(0,n,numDim)] - X0;
114 register const double D1=mesh-> Nodes->Coordinates[INDEX2(1,n,numDim)] - X1;
115 register const double d= D0 * D0 + D1 * D1;
116 if ( d < dist_local) {
117 dist_local = d;
118 node_id_local=n;
119 }
120 }
121 #pragma omp critical
122 {
123 if ( dist_local < dist_p[i] ) {
124 dist_p[i] = dist_local;
125 node_id_p[i] = node_id_local;
126 }
127 }
128 }
129 }
130 } else {
131 #pragma omp parallel private(i)
132 {
133 for (i=0; i< numPoints; ++i) {
134 register const double X0=points_ptr[INDEX2(0,i,numDim)];
135 register double dist_local=LARGE_POSITIVE_FLOAT;
136 register index_t node_id_local=-1;
137 #pragma omp for private(n)
138 for (n=0; n< mesh-> Nodes->numNodes; n++) {
139 register const double D0=mesh-> Nodes->Coordinates[INDEX2(0,n,numDim)] - X0;
140 register const double d= D0 * D0;
141 if ( d < dist_local) {
142 dist_local = d;
143 node_id_local=n;
144 }
145 }
146 #pragma omp critical
147 {
148 if ( dist_local < dist_p[i] ) {
149 dist_p[i] = dist_local;
150 node_id_p[i] = node_id_local;
151 }
152 }
153 }
154 }
155 }
156
157 /* now we need to reduce this across all processors */
158 #ifdef ESYS_MPI
159 {
160 MPI_Op op;
161 int count = 2* numPoints;
162 double *sendbuf =NULL, *recvbuf=NULL;
163 sendbuf=new double[count];
164 recvbuf=new double[count];
165
166 for (i=0; i< numPoints; ++i) {
167 sendbuf[2*i ]= dist_p[i];
168 sendbuf[2*i+1]= (double) (mesh-> Nodes->Id[node_id_p[i]]);
169 }
170 MPI_Op_create(Finley_Mesh_MPI_minimizeDistance, TRUE,&op);
171 MPI_Allreduce ( sendbuf, recvbuf, count, MPI_DOUBLE, op, mpi_info->comm);
172 MPI_Op_free(&op);
173 /* if the node id has changed we found another node which is closer elsewhere */
174 for (i=0; i< numPoints; ++i) {
175 register const int best_fit_Id = (int) (recvbuf[2*i+1] + 0.5);
176 if ( best_fit_Id != mesh-> Nodes->Id[node_id_p[i]] ) {
177 node_id_p[i] =-1;
178 }
179 }
180 delete[] sendbuf;
181 delete[] recvbuf;
182 }
183 #endif
184 /* we pick the points to be used on this processor */
185 numNewPoints=0;
186 for (i=0; i< numPoints; ++i) {
187 if (node_id_p[i]>-1) { /* this processor uses a node which is identical to point i */
188
189 if ( mesh-> Nodes->globalReducedDOFIndex[node_id_p[i]] >-1) { /* the point is also used in the reduced mesh */
190
191 register const index_t global_id=mesh->Nodes->globalDegreesOfFreedom[node_id_p[i]];
192 if ( (firstDOF<= global_id) && ( global_id <lastDOF) ) { /* is this point actually relevant */
193
194
195 /* is this point already in the Point table? */
196 bool_t notAnOldPoint=TRUE;
197
198 if (numOldPoints >0) {
199 for (k=0; k<numOldPoints; ++k) {
200 if ( global_id == oldPoints->Nodes[k]) {
201 notAnOldPoint=FALSE;
202 break;
203 }
204 }
205 }
206 if (notAnOldPoint) {
207 /* is this point unique in the new list of points? */
208 bool_t notANewPoint=TRUE;
209 for (k=0; k<numNewPoints; ++k) {
210 if ( global_id == mesh->Nodes->globalDegreesOfFreedom[node_id_p[point_index_p[k]]]) {
211 notANewPoint=FALSE;
212 break;
213 }
214 }
215 if (notANewPoint) {
216 point_index_p[numNewPoints]=i;
217 numNewPoints++;
218 }
219 }
220 }
221 }
222 }
223 }
224 /* now we are ready to create the new Point table */
225 Finley_ElementFile_allocTable(newPoints,numOldPoints+numNewPoints);
226 if (numOldPoints > 0) {
227 #pragma omp parallel for private(n) schedule(static)
228 for(n=0;n<numOldPoints;n++) {
229 newPoints->Owner[n]=oldPoints->Owner[n];
230 newPoints->Id[n] =oldPoints->Id[n];
231 newPoints->Tag[n] =oldPoints->Tag[n];
232 newPoints->Nodes[n]=oldPoints->Nodes[n];
233 newPoints->Color[n]=0;
234 }
235 }
236 #pragma omp parallel for private(n) schedule(static)
237 for(n=0;n<numNewPoints;n++) {
238 register const index_t idx = point_index_p[n];
239 newPoints->Owner[numOldPoints+n]=mpi_info->rank;
240 newPoints->Id[numOldPoints+n] =0;
241 newPoints->Tag[numOldPoints+n] =tags_ptr[idx];
242 newPoints->Nodes[numOldPoints+n]=node_id_p[idx];
243 newPoints->Color[numOldPoints+n]=0;
244 }
245 newPoints->minColor=0;
246 newPoints->maxColor=0;
247
248 /* all done */
249 delete[] dist_p;
250 delete[] node_id_p;
251 delete[] point_index_p;
252 Finley_ReferenceElementSet_dealloc(refPoints);
253 Esys_MPIInfo_free(mpi_info);
254 if (Finley_noError()) {
255 Finley_ElementFile_free(oldPoints);
256 mesh->Points=newPoints;
257 } else {
258 Finley_ElementFile_free(newPoints);
259 }
260 }

  ViewVC Help
Powered by ViewVC 1.1.26