/[escript]/trunk/dudley/src/Assemble_getSize.c
ViewVC logotype

Contents of /trunk/dudley/src/Assemble_getSize.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File MIME type: text/plain
File size: 3755 byte(s)
Merging dudley and scons updates from branches

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 /* assemblage routines: */
17
18 /* calculates the minimum distance between two vertices of elements and assigns the value to each */
19 /* quadrature point in element_size */
20
21 /**************************************************************/
22
23 #include "Assemble.h"
24 #include "Util.h"
25 #ifdef _OPENMP
26 #include <omp.h>
27 #endif
28
29 /**************************************************************/
30 void Dudley_Assemble_getSize(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, escriptDataC * element_size)
31 {
32
33 double *local_X = NULL, *element_size_array;
34 dim_t e, n0, n1, q, i, NVertices, NN, NS, numQuad, numDim;
35 double d, diff, min_diff;
36 Dudley_resetError();
37
38 if (nodes == NULL || elements == NULL)
39 {
40 return;
41 }
42
43 numDim = nodes->numDim;
44
45 /* now we look up what type of elements we need based on the function space of element_size */
46 /* if it is DUDLEY_REDUCED_ELEMENTS or DUDLEY_REDUCED_FACE_ELEMENTS then we have single quad point */
47
48 if (Dudley_Assemble_reducedIntegrationOrder(element_size))
49 {
50 numQuad = 1;
51 }
52 else
53 {
54 numQuad = elements->numDim + 1;
55 }
56
57 NN = elements->numNodes;
58 NS = elements->numDim + 1;
59 NVertices = elements->numDim + 1;
60
61 /* check the dimensions of element_size */
62
63 if (!numSamplesEqual(element_size, numQuad, elements->numElements))
64 {
65 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getSize: illegal number of samples of element size Data object");
66 }
67 else if (!isDataPointShapeEqual(element_size, 0, &(numDim)))
68 {
69 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getSize: illegal data point shape of element size Data object");
70 }
71 else if (!isExpanded(element_size))
72 {
73 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getSize: expanded Data object is expected for element size.");
74 }
75 /* now we can start: */
76
77 if (Dudley_noError())
78 {
79 requireWrite(element_size);
80 #pragma omp parallel private(local_X)
81 {
82 /* allocation of work arrays */
83 local_X = THREAD_MEMALLOC(NN * numDim, double);
84 if (!Dudley_checkPtr(local_X))
85 {
86 /* open the element loop */
87 #pragma omp for private(e,min_diff,diff,n0,n1,d,q,i,element_size_array) schedule(static)
88 for (e = 0; e < elements->numElements; e++)
89 {
90 /* gather local coordinates of nodes into local_X(numDim,NN): */
91 Dudley_Util_Gather_double(NS, &(elements->Nodes[INDEX2(0, e, NN)]), numDim, nodes->Coordinates,
92 local_X);
93 /* calculate minimal differences */
94 min_diff = -1;
95 for (n0 = 0; n0 < NVertices; n0++)
96 {
97 for (n1 = n0 + 1; n1 < NVertices; n1++)
98 {
99 diff = 0;
100 for (i = 0; i < numDim; i++)
101 {
102 d = local_X[INDEX2(i, n0, numDim)] - local_X[INDEX2(i, n1, numDim)];
103 diff += d * d;
104 }
105 if (min_diff < 0)
106 {
107 min_diff = diff;
108 }
109 else
110 {
111 min_diff = MIN(min_diff, diff);
112 }
113 }
114 }
115 min_diff = sqrt(MAX(min_diff, 0));
116 /* set all values to min_diff */
117 element_size_array = getSampleDataRW(element_size, e);
118 for (q = 0; q < numQuad; q++)
119 element_size_array[q] = min_diff;
120 }
121 }
122 THREAD_MEMFREE(local_X);
123 } /* end of parallel region */
124 }
125 return;
126 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26