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

Contents of /trunk/finley/src/Assemble_interpolate.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: 7050 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: interpolates nodal data in a data array onto elements
22 (=integration points)
23
24 *****************************************************************************/
25
26 #include "Assemble.h"
27 #include "Util.h"
28
29 #include <escript/index.h>
30
31 namespace finley {
32
33 template<typename Scalar>
34 void Assemble_interpolate(const NodeFile* nodes, const ElementFile* elements,
35 const escript::Data& data,
36 escript::Data& interpolated_data)
37 {
38 if (!nodes || !elements)
39 return;
40
41 const int data_type = data.getFunctionSpace().getTypeCode();
42 const bool reducedOrder = util::hasReducedIntegrationOrder(interpolated_data);
43 const_ReferenceElement_ptr refElement(elements->referenceElementSet->
44 borrowReferenceElement(reducedOrder));
45
46 const int* resort_nodes = NULL;
47 const index_t* map = NULL;
48 int numSub = 0;
49 dim_t numNodes = 0;
50 const_ShapeFunction_ptr basis;
51 int dof_offset = 0;
52
53 if (data_type == FINLEY_NODES) {
54 numSub = refElement->Type->numSubElements;
55 resort_nodes=refElement->Type->subElementNodes;
56 basis=refElement->BasisFunctions;
57 numNodes=nodes->getNumNodes();
58 map=nodes->borrowTargetNodes();
59 if (interpolated_data.getFunctionSpace().getTypeCode()==FINLEY_CONTACT_ELEMENTS_2) {
60 dof_offset=refElement->Type->offsets[1];
61 } else {
62 dof_offset=refElement->Type->offsets[0];
63 }
64 } else if (data_type==FINLEY_REDUCED_NODES) {
65 numSub=1;
66 resort_nodes=refElement->Type->linearNodes;
67 basis=refElement->LinearBasisFunctions;
68 numNodes=nodes->getNumReducedNodes();
69 map=nodes->borrowTargetReducedNodes();
70 if (interpolated_data.getFunctionSpace().getTypeCode()==FINLEY_CONTACT_ELEMENTS_2) {
71 dof_offset=refElement->LinearType->offsets[1];
72 } else {
73 dof_offset=refElement->LinearType->offsets[0];
74 }
75 } else if (data_type == FINLEY_DEGREES_OF_FREEDOM) {
76 if (elements->MPIInfo->size > 1) {
77 throw escript::ValueError("Assemble_interpolate: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
78 }
79 numSub = refElement->Type->numSubElements;
80 resort_nodes = refElement->Type->subElementNodes;
81 basis = refElement->BasisFunctions;
82 numNodes = nodes->getNumDegreesOfFreedom();
83 map = nodes->borrowTargetDegreesOfFreedom();
84 if (interpolated_data.getFunctionSpace().getTypeCode()==FINLEY_CONTACT_ELEMENTS_2) {
85 dof_offset = refElement->Type->offsets[1];
86 } else {
87 dof_offset = refElement->Type->offsets[0];
88 }
89 } else if (data_type == FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
90 if (elements->MPIInfo->size > 1) {
91 throw escript::ValueError("Assemble_interpolate: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
92 }
93 numSub=1;
94 resort_nodes=refElement->Type->linearNodes;
95 basis=refElement->LinearBasisFunctions;
96 numNodes=nodes->getNumReducedDegreesOfFreedom();
97 map=nodes->borrowTargetReducedDegreesOfFreedom();
98 if (interpolated_data.getFunctionSpace().getTypeCode()==FINLEY_CONTACT_ELEMENTS_2) {
99 dof_offset=refElement->LinearType->offsets[1];
100 } else {
101 dof_offset=refElement->LinearType->offsets[0];
102 }
103 } else {
104 throw escript::ValueError("Assemble_interpolate: invalid functionspace");
105 }
106
107 const int numComps = data.getDataPointSize();
108 const int numQuad = basis->numQuadNodes;
109 const int numShapesTotal = basis->Type->numShapes*refElement->Type->numSides;
110 const int NN = elements->numNodes;
111 const int NS_DOF = basis->Type->numShapes;
112
113 // check the dimensions of interpolated_data and data
114 if (!interpolated_data.numSamplesEqual(numQuad*numSub, elements->numElements)) {
115 throw escript::ValueError("Assemble_interpolate: illegal number of samples of output Data object");
116 } else if (!data.numSamplesEqual(1,numNodes)) {
117 throw escript::ValueError("Assemble_interpolate: illegal number of samples of input Data object");
118 } else if (numComps != interpolated_data.getDataPointSize()) {
119 throw escript::ValueError("Assemble_interpolate: number of components of input and interpolated Data do not match.");
120 } else if (!interpolated_data.actsExpanded()) {
121 throw escript::ValueError("Assemble_interpolate: expanded Data object is expected for output data.");
122 }
123
124 const Scalar zero = static_cast<Scalar>(0);
125 interpolated_data.requireWrite();
126 #pragma omp parallel
127 {
128 std::vector<Scalar> local_data(NS_DOF * numComps * numSub);
129 const size_t numComps_size = numComps * sizeof(Scalar);
130 // open the element loop
131 #pragma omp for
132 for (index_t e = 0; e < elements->numElements; e++) {
133 for (int isub = 0; isub < numSub; isub++) {
134 for (int q = 0; q < NS_DOF; q++) {
135 const index_t i = elements->Nodes[INDEX2(resort_nodes[INDEX2(dof_offset+q,isub,numShapesTotal)],e,NN)];
136 const Scalar* data_array = data.getSampleDataRO(map[i], zero);
137 memcpy(&local_data[INDEX3(0, q, isub, numComps,NS_DOF)], data_array, numComps_size);
138 }
139 }
140 // calculate interpolated_data=local_data*S
141 util::smallMatSetMult1<Scalar>(numSub, numComps, numQuad,
142 interpolated_data.getSampleDataRW(e, zero), NS_DOF,
143 local_data, basis->S);
144 } // end of element loop
145 } // end of parallel region
146 }
147
148 // instantiate our two supported versions
149 template void Assemble_interpolate<escript::DataTypes::real_t>(
150 const NodeFile* nodes, const ElementFile* elements,
151 const escript::Data& data,
152 escript::Data& interpolated_data);
153 template void Assemble_interpolate<escript::DataTypes::cplx_t>(
154 const NodeFile* nodes, const ElementFile* elements,
155 const escript::Data& data,
156 escript::Data& interpolated_data);
157
158 } // namespace finley
159

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26