/[escript]/trunk/finley/src/Assemble_AverageElementData.cpp
ViewVC logotype

Contents of /trunk/finley/src/Assemble_AverageElementData.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 4790 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18
19 /****************************************************************************
20
21 Assemblage routines: averages data between elements.
22
23 *****************************************************************************/
24
25 #include "Assemble.h"
26 #include "Util.h"
27
28 #include <escript/index.h>
29
30 namespace finley {
31
32 template<typename Scalar>
33 void Assemble_AverageElementData(const ElementFile* elements,
34 escript::Data& out, const escript::Data& in)
35 {
36 if (!elements)
37 return;
38
39 const double* wq;
40 int numQuad_in, numQuad_out;
41 if (util::hasReducedIntegrationOrder(in)) {
42 numQuad_in = elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
43 wq = &elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->QuadWeights[0];
44 } else {
45 numQuad_in = elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
46 wq = &elements->referenceElementSet->referenceElement->Parametrization->QuadWeights[0];
47 }
48 if (util::hasReducedIntegrationOrder(out)) {
49 numQuad_out = elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
50 } else {
51 numQuad_out = elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
52 }
53
54 // check out and in
55 const dim_t numElements = elements->numElements;
56 const int numComps = out.getDataPointSize();
57 const Scalar zero = static_cast<Scalar>(0);
58
59 if (numComps != in.getDataPointSize()) {
60 throw escript::ValueError("Assemble_AverageElementData: number of components of input and output data do not match.");
61 } else if (!in.numSamplesEqual(numQuad_in,numElements)) {
62 throw escript::ValueError("Assemble_AverageElementData: illegal number of samples of input Data object");
63 } else if (!out.numSamplesEqual(numQuad_out,numElements)) {
64 throw escript::ValueError("Assemble_AverageElementData: illegal number of samples of output Data object");
65 } else if (!out.actsExpanded()) {
66 throw escript::ValueError("Assemble_AverageElementData: expanded Data object is expected for output data.");
67 } else if (in.isComplex() != out.isComplex()) {
68 throw escript::ValueError("Assemble_AverageElementData: complexity of input and output data must match.");
69 } else {
70 out.requireWrite();
71 if (in.actsExpanded()) {
72 double vol = 0.;
73 for (int q = 0; q < numQuad_in; ++q)
74 vol += wq[q];
75 const double volinv = 1./vol;
76 #pragma omp parallel for
77 for (index_t n = 0; n < numElements; n++) {
78 const Scalar* in_array = in.getSampleDataRO(n, zero);
79 Scalar* out_array = out.getSampleDataRW(n, zero);
80 for (int i = 0; i < numComps; ++i) {
81 Scalar rtmp = zero;
82 for (int q = 0; q < numQuad_in; ++q)
83 rtmp += in_array[INDEX2(i,q,numComps)] * wq[q];
84 rtmp *= volinv;
85 for (int q = 0; q < numQuad_out; ++q)
86 out_array[INDEX2(i,q,numComps)] = rtmp;
87 }
88 }
89 } else { // constant data
90 const size_t numComps_size = numComps * sizeof(Scalar);
91 #pragma omp parallel for
92 for (index_t n = 0; n < numElements; n++) {
93 const Scalar* in_array = in.getSampleDataRO(n, zero);
94 Scalar* out_array = out.getSampleDataRW(n, zero);
95 for (int q = 0; q < numQuad_out; q++)
96 memcpy(out_array+q*numComps, in_array, numComps_size);
97 }
98 }
99 }
100 }
101
102 // instantiate our two supported versions
103 template void Assemble_AverageElementData<escript::DataTypes::real_t>(
104 const ElementFile* elements,
105 escript::Data& out, const escript::Data& in);
106 template void Assemble_AverageElementData<escript::DataTypes::cplx_t>(
107 const ElementFile* elements,
108 escript::Data& out, const escript::Data& in);
109
110 } // namespace finley
111

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/finley/src/Assemble_AverageElementData.cpp:5567-5588 /branches/lapack2681/finley/src/Assemble_AverageElementData.cpp:2682-2741 /branches/pasowrap/finley/src/Assemble_AverageElementData.cpp:3661-3674 /branches/py3_attempt2/finley/src/Assemble_AverageElementData.cpp:3871-3891 /branches/restext/finley/src/Assemble_AverageElementData.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Assemble_AverageElementData.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_AverageElementData.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Assemble_AverageElementData.cpp:3471-3974 /branches/trilinos_from_5897/finley/src/Assemble_AverageElementData.cpp:5898-6118 /release/3.0/finley/src/Assemble_AverageElementData.cpp:2591-2601 /release/4.0/finley/src/Assemble_AverageElementData.cpp:5380-5406 /trunk/finley/src/Assemble_AverageElementData.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26