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

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26