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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4154 - (show annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 3784 byte(s)
Round 1 of copyright fixes
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: copies data between elements */
19
20 /************************************************************************************/
21
22 #include "Assemble.h"
23 #include "Util.h"
24 #ifdef _OPENMP
25 #include <omp.h>
26 #endif
27 /****************************************************************************************************************************/
28
29 #include "ShapeTable.h"
30
31 void Dudley_Assemble_AverageElementData(Dudley_ElementFile * elements, escriptDataC * out, escriptDataC * in)
32 {
33 dim_t n, q, numElements, numQuad_in, numQuad_out, i;
34 __const double *in_array;
35 double *out_array, vol, volinv, wq;
36 register double rtmp;
37 dim_t numComps = getDataPointSize(out);
38 size_t numComps_size;
39
40 Dudley_resetError();
41 if (elements == NULL)
42 {
43 return;
44 }
45
46 numElements = elements->numElements;
47 if (Dudley_Assemble_reducedIntegrationOrder(in))
48 {
49 numQuad_in = QuadNums[elements->numDim][0];
50 wq = QuadWeight[elements->numDim][0];
51
52 }
53 else
54 {
55 numQuad_in = QuadNums[elements->numDim][1];
56 wq = QuadWeight[elements->numDim][1];
57 }
58 if (Dudley_Assemble_reducedIntegrationOrder(out))
59 {
60 numQuad_out = QuadNums[elements->numDim][0];
61 }
62 else
63 {
64 numQuad_out = QuadNums[elements->numDim][1];
65
66 }
67
68 /* check out and in */
69 if (numComps != getDataPointSize(in))
70 {
71 Dudley_setError(TYPE_ERROR,
72 "Dudley_Assemble_AverageElementData: number of components of input and output Data do not match.");
73 }
74 else if (!numSamplesEqual(in, numQuad_in, numElements))
75 {
76 Dudley_setError(TYPE_ERROR,
77 "Dudley_Assemble_AverageElementData: illegal number of samples of input Data object");
78 }
79 else if (!numSamplesEqual(out, numQuad_out, numElements))
80 {
81 Dudley_setError(TYPE_ERROR,
82 "Dudley_Assemble_AverageElementData: illegal number of samples of output Data object");
83 }
84 else if (!isExpanded(out))
85 {
86 Dudley_setError(TYPE_ERROR,
87 "Dudley_Assemble_AverageElementData: expanded Data object is expected for output data.");
88 }
89
90 /* now we can start */
91
92 if (Dudley_noError())
93 {
94 if (isExpanded(in))
95 {
96 vol = 0;
97 for (q = 0; q < numQuad_in; ++q)
98 vol += wq;
99 volinv = 1. / vol;
100 requireWrite(out);
101 #pragma omp parallel private(n, i, rtmp, q, in_array, out_array)
102 {
103 # pragma omp for schedule(static)
104 for (n = 0; n < numElements; n++)
105 {
106 in_array = getSampleDataRO(in, n);
107 out_array = getSampleDataRW(out, n);
108 for (i = 0; i < numComps; ++i)
109 {
110 rtmp = 0;
111 for (q = 0; q < numQuad_in; ++q)
112 rtmp += in_array[INDEX2(i, q, numComps)] * wq;
113 rtmp *= volinv;
114 for (q = 0; q < numQuad_out; ++q)
115 out_array[INDEX2(i, q, numComps)] = rtmp;
116 }
117 }
118 }
119 }
120 else
121 {
122 numComps_size = numComps * sizeof(double);
123 requireWrite(out);
124 #pragma omp parallel private(q,n,out_array,in_array)
125 {
126 # pragma omp for schedule(static)
127 for (n = 0; n < numElements; n++)
128 {
129 in_array = getSampleDataRO(in, n);
130 out_array = getSampleDataRW(out, n);
131 for (q = 0; q < numQuad_out; q++)
132 memcpy(out_array + q * numComps, in_array, numComps_size);
133 }
134 }
135 }
136 }
137 return;
138 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26