/[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 4492 - (show annotations)
Tue Jul 2 01:44:11 2013 UTC (5 years, 9 months ago) by caltinay
File size: 6589 byte(s)
Finley changes that were held back while in release mode - moved more stuff
into finley namespace.

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /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 /release/3.0/finley/src/Assemble_interpolate.cpp:2591-2601 /trunk/finley/src/Assemble_interpolate.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26