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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years, 3 months ago) by jfenwick
File MIME type: text/plain
File size: 3191 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: integrates data on quadrature points */
17
18 /**************************************************************/
19
20 #include "Assemble.h"
21 #include "Util.h"
22 #ifdef _OPENMP
23 #include <omp.h>
24 #endif
25
26 /**************************************************************/
27
28 void Dudley_Assemble_integrate(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, escriptDataC * data, double *out)
29 {
30 /* type_t data_type=getFunctionSpaceType(data);*/
31 dim_t numQuadTotal;
32 dim_t numComps = getDataPointSize(data);
33 Dudley_ElementFile_Jacobeans *jac = NULL;
34 Esys_MPI_rank my_mpi_rank;
35
36 Dudley_resetError();
37 if (nodes == NULL || elements == NULL)
38 return;
39 my_mpi_rank = nodes->MPIInfo->rank;
40 /* set some parameter */
41 jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, Dudley_Assemble_reducedIntegrationOrder(data));
42 if (Dudley_noError())
43 {
44 numQuadTotal = jac->numQuad;
45 /* check the shape of the data */
46 if (!numSamplesEqual(data, numQuadTotal, elements->numElements))
47 {
48 Dudley_setError(TYPE_ERROR,
49 "Dudley_Assemble_integrate: illegal number of samples of integrant kernel Data object");
50 }
51 /* now we can start */
52
53 if (Dudley_noError())
54 {
55 dim_t q, e, i;
56 __const double *data_array = NULL;
57 double *out_local = NULL, rtmp;
58 for (q = 0; q < numComps; q++)
59 out[q] = 0;
60 #pragma omp parallel private(q,i,rtmp,data_array,out_local)
61 {
62 out_local = THREAD_MEMALLOC(numComps, double);
63 if (!Dudley_checkPtr(out_local))
64 {
65 /* initialize local result */
66
67 for (i = 0; i < numComps; i++)
68 out_local[i] = 0;
69
70 /* open the element loop */
71
72 if (isExpanded(data))
73 {
74 #pragma omp for private(e) schedule(static)
75 for (e = 0; e < elements->numElements; e++)
76 {
77 if (elements->Owner[e] == my_mpi_rank)
78 {
79 double vol = jac->absD[e] * jac->quadweight;
80 data_array = getSampleDataRO(data, e);
81 for (q = 0; q < numQuadTotal; q++)
82 {
83 for (i = 0; i < numComps; i++)
84 out_local[i] += data_array[INDEX2(i, q, numComps)] * vol;
85 }
86 }
87 }
88 }
89 else
90 {
91 #pragma omp for private(e) schedule(static)
92 for (e = 0; e < elements->numElements; e++)
93 {
94 if (elements->Owner[e] == my_mpi_rank)
95 {
96 double vol = jac->absD[e] * jac->quadweight;
97 data_array = getSampleDataRO(data, e);
98 rtmp = 0.;
99 for (q = 0; q < numQuadTotal; q++)
100 rtmp += vol;
101 for (i = 0; i < numComps; i++)
102 out_local[i] += data_array[i] * rtmp;
103 }
104 }
105 }
106 /* add local results to global result */
107 #pragma omp critical
108 for (i = 0; i < numComps; i++)
109 out[i] += out_local[i];
110 }
111 THREAD_MEMFREE(out_local);
112 }
113 }
114 }
115 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26