/[escript]/branches/doubleplusgood/dudley/src/Assemble_getSize.cpp
ViewVC logotype

Annotation of /branches/doubleplusgood/dudley/src/Assemble_getSize.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4332 - (hide annotations)
Thu Mar 21 04:21:14 2013 UTC (6 years ago) by jfenwick
File size: 3833 byte(s)
like sand though the hourglass
1 jgs 82
2 jfenwick 3981 /*****************************************************************************
3 ksteube 1811 *
4 jfenwick 4154 * Copyright (c) 2003-2013 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 ksteube 1811 *
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 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     * Development since 2012 by School of Earth Sciences
13     *
14     *****************************************************************************/
15 ksteube 1312
16 jfenwick 3981 /************************************************************************************/
17 jgs 82
18     /* assemblage routines: */
19    
20     /* calculates the minimum distance between two vertices of elements and assigns the value to each */
21     /* quadrature point in element_size */
22    
23 jfenwick 3981 /************************************************************************************/
24 jgs 82
25     #include "Assemble.h"
26     #include "Util.h"
27     #ifdef _OPENMP
28     #include <omp.h>
29     #endif
30    
31 jfenwick 3981 /************************************************************************************/
32 jfenwick 3224 void Dudley_Assemble_getSize(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, escriptDataC * element_size)
33 jfenwick 3199 {
34 jgs 82
35 jfenwick 3224 double *local_X = NULL, *element_size_array;
36     dim_t e, n0, n1, q, i, NVertices, NN, NS, numQuad, numDim;
37 gross 3805 double d, diff, max_diff;
38 jfenwick 3224 Dudley_resetError();
39 jgs 82
40 jfenwick 3224 if (nodes == NULL || elements == NULL)
41     {
42 jfenwick 3199 return;
43 jfenwick 3224 }
44 jfenwick 3199
45 jfenwick 3224 numDim = nodes->numDim;
46 jfenwick 3199
47 jfenwick 3251 /* now we look up what type of elements we need based on the function space of element_size */
48     /* if it is DUDLEY_REDUCED_ELEMENTS or DUDLEY_REDUCED_FACE_ELEMENTS then we have single quad point */
49 jfenwick 3199
50 jfenwick 3224 if (Dudley_Assemble_reducedIntegrationOrder(element_size))
51     {
52     numQuad = 1;
53     }
54     else
55     {
56     numQuad = elements->numDim + 1;
57     }
58 jfenwick 3199
59 jfenwick 3224 NN = elements->numNodes;
60     NS = elements->numDim + 1;
61     NVertices = elements->numDim + 1;
62 jgs 82
63 jfenwick 3224 /* check the dimensions of element_size */
64 jgs 82
65 jfenwick 3224 if (!numSamplesEqual(element_size, numQuad, elements->numElements))
66     {
67     Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getSize: illegal number of samples of element size Data object");
68     }
69     else if (!isDataPointShapeEqual(element_size, 0, &(numDim)))
70     {
71     Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getSize: illegal data point shape of element size Data object");
72     }
73     else if (!isExpanded(element_size))
74     {
75     Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getSize: expanded Data object is expected for element size.");
76     }
77     /* now we can start: */
78    
79     if (Dudley_noError())
80     {
81     requireWrite(element_size);
82     #pragma omp parallel private(local_X)
83     {
84     /* allocation of work arrays */
85 jfenwick 4332 local_X = new double[NN * numDim];
86 jfenwick 3224 if (!Dudley_checkPtr(local_X))
87     {
88     /* open the element loop */
89 gross 3805 #pragma omp for private(e,max_diff,diff,n0,n1,d,q,i,element_size_array) schedule(static)
90 jfenwick 3224 for (e = 0; e < elements->numElements; e++)
91     {
92     /* gather local coordinates of nodes into local_X(numDim,NN): */
93     Dudley_Util_Gather_double(NS, &(elements->Nodes[INDEX2(0, e, NN)]), numDim, nodes->Coordinates,
94     local_X);
95     /* calculate minimal differences */
96 gross 3805 max_diff = 0;
97 jfenwick 3224 for (n0 = 0; n0 < NVertices; n0++)
98     {
99     for (n1 = n0 + 1; n1 < NVertices; n1++)
100     {
101     diff = 0;
102     for (i = 0; i < numDim; i++)
103     {
104     d = local_X[INDEX2(i, n0, numDim)] - local_X[INDEX2(i, n1, numDim)];
105     diff += d * d;
106     }
107 gross 3805
108     max_diff = MAX(max_diff, diff);
109    
110 gross 2748 }
111 jfenwick 3224 }
112 gross 3805 max_diff = sqrt(max_diff);
113     /* set all values to max_diff */
114 jfenwick 3224 element_size_array = getSampleDataRW(element_size, e);
115     for (q = 0; q < numQuad; q++)
116 gross 3805 element_size_array[q] = max_diff;
117 jfenwick 3224 }
118     }
119 jfenwick 4332 delete[] local_X;
120 jfenwick 3224 } /* end of parallel region */
121     }
122     return;
123 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26