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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4332 - (show annotations)
Thu Mar 21 04:21:14 2013 UTC (5 years, 11 months ago) by jfenwick
File size: 3833 byte(s)
like sand though the hourglass
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 /* 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 /************************************************************************************/
24
25 #include "Assemble.h"
26 #include "Util.h"
27 #ifdef _OPENMP
28 #include <omp.h>
29 #endif
30
31 /************************************************************************************/
32 void Dudley_Assemble_getSize(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, escriptDataC * element_size)
33 {
34
35 double *local_X = NULL, *element_size_array;
36 dim_t e, n0, n1, q, i, NVertices, NN, NS, numQuad, numDim;
37 double d, diff, max_diff;
38 Dudley_resetError();
39
40 if (nodes == NULL || elements == NULL)
41 {
42 return;
43 }
44
45 numDim = nodes->numDim;
46
47 /* 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
50 if (Dudley_Assemble_reducedIntegrationOrder(element_size))
51 {
52 numQuad = 1;
53 }
54 else
55 {
56 numQuad = elements->numDim + 1;
57 }
58
59 NN = elements->numNodes;
60 NS = elements->numDim + 1;
61 NVertices = elements->numDim + 1;
62
63 /* check the dimensions of element_size */
64
65 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 local_X = new double[NN * numDim];
86 if (!Dudley_checkPtr(local_X))
87 {
88 /* open the element loop */
89 #pragma omp for private(e,max_diff,diff,n0,n1,d,q,i,element_size_array) schedule(static)
90 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 max_diff = 0;
97 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
108 max_diff = MAX(max_diff, diff);
109
110 }
111 }
112 max_diff = sqrt(max_diff);
113 /* set all values to max_diff */
114 element_size_array = getSampleDataRW(element_size, e);
115 for (q = 0; q < numQuad; q++)
116 element_size_array[q] = max_diff;
117 }
118 }
119 delete[] local_X;
120 } /* end of parallel region */
121 }
122 return;
123 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26