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

Contents of /trunk/finley/src/Assemble_PDE_System_C.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: 9229 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 Assembles the system of numEqu PDEs into the stiffness matrix S and right
21 hand side F
22
23 D_{k,m} u_m and Y_k
24
25 u has p.numComp components in a 3D domain. The shape functions for test and
26 solution must be identical and 2*row_NS == row_NN (contact elements).
27
28 Shape of the coefficients:
29
30 D = p.numEqu x p.numComp
31 Y = p.numEqu
32
33 *****************************************************************************/
34
35 #include "Assemble.h"
36 #include "Util.h"
37
38 #include <escript/index.h>
39
40 namespace finley {
41
42 template<typename Scalar>
43 void Assemble_PDE_System_C(const AssembleParameters& p, const escript::Data& D,
44 const escript::Data& Y)
45 {
46 bool expandedD = D.actsExpanded();
47 bool expandedY = Y.actsExpanded();
48 const Scalar zero = static_cast<Scalar>(0);
49 Scalar* F_p = NULL;
50 if (!p.F.isEmpty()) {
51 p.F.requireWrite();
52 F_p = p.F.getSampleDataRW(0, zero);
53 }
54 const std::vector<double>& S(p.row_jac->BasisFunctions->S);
55
56 #pragma omp parallel
57 {
58 std::vector<Scalar> EM_S(p.row_numShapesTotal*p.col_numShapesTotal*p.numEqu*p.numComp);
59 std::vector<Scalar> EM_F(p.row_numShapesTotal*p.numEqu);
60 IndexVector row_index(p.row_numShapesTotal);
61
62 for (index_t color = p.elements->minColor; color <= p.elements->maxColor; color++) {
63 // loop over all elements
64 #pragma omp for
65 for (index_t e = 0; e < p.elements->numElements; e++) {
66 if (p.elements->Color[e] == color) {
67 for (int isub = 0; isub < p.numSub; isub++) {
68 const double* Vol = &(p.row_jac->volume[INDEX3(0,isub,e,p.numQuadSub,p.numSub)]);
69 bool add_EM_F = false;
70 bool add_EM_S = false;
71 ///////////////
72 // process D //
73 ///////////////
74 if (!D.isEmpty()) {
75 const Scalar* D_p = D.getSampleDataRO(e, zero);
76 add_EM_S = true;
77 if (expandedD) {
78 const Scalar* D_q = &D_p[INDEX4(0,0,0,isub,p.numEqu,p.numComp, p.numQuadSub)];
79 for (int s = 0; s < p.row_numShapes; s++) {
80 for (int r = 0; r < p.col_numShapes; r++) {
81 for (int k = 0; k < p.numEqu; k++) {
82 for (int m = 0; m < p.numComp; m++) {
83 Scalar val = zero;
84 for (int q = 0; q < p.numQuadSub; q++) {
85 val += Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_q[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_numShapes)];
86 }
87 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]= val;
88 EM_S[INDEX4(k,m,s,r+p.col_numShapes,p.numEqu,p.numComp,p.row_numShapesTotal)]=-val;
89 EM_S[INDEX4(k,m,s+p.row_numShapes,r,p.numEqu,p.numComp,p.row_numShapesTotal)]=-val;
90 EM_S[INDEX4(k,m,s+p.row_numShapes,r+p.col_numShapes,p.numEqu,p.numComp,p.row_numShapesTotal)]= val;
91 }
92 }
93 }
94 }
95 } else { // constant D
96 for (int s = 0; s < p.row_numShapes; s++) {
97 for (int r = 0; r < p.col_numShapes; r++) {
98 Scalar f = zero;
99 for (int q = 0; q < p.numQuadSub; q++) {
100 f += Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];
101 }
102 for (int k = 0; k < p.numEqu; k++) {
103 for (int m = 0; m < p.numComp; m++) {
104 const Scalar fD = f * D_p[INDEX2(k,m,p.numEqu)];
105 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]= fD;
106 EM_S[INDEX4(k,m,s,r+p.col_numShapes,p.numEqu,p.numComp,p.row_numShapesTotal)]=-fD;
107 EM_S[INDEX4(k,m,s+p.row_numShapes,r,p.numEqu,p.numComp,p.row_numShapesTotal)]=-fD;
108 EM_S[INDEX4(k,m,s+p.row_numShapes,r+p.col_numShapes,p.numEqu,p.numComp,p.row_numShapesTotal)]= fD;
109 }
110 }
111 }
112 }
113 }
114 }
115 ///////////////
116 // process Y //
117 ///////////////
118 if (!Y.isEmpty()) {
119 const Scalar* Y_p = Y.getSampleDataRO(e, zero);
120 add_EM_F = true;
121 if (expandedY) {
122 const Scalar* Y_q = &Y_p[INDEX3(0, 0, isub, p.numEqu, p.numQuadSub)];
123 for (int s = 0; s < p.row_numShapes; s++) {
124 for (int k = 0; k < p.numEqu; k++) {
125 Scalar val = zero;
126 for (int q = 0; q < p.numQuadSub; q++)
127 val += Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*Y_q[INDEX2(k,q,p.numEqu)];
128 EM_F[INDEX2(k,s,p.numEqu)] = -val;
129 EM_F[INDEX2(k,s+p.row_numShapes,p.numEqu)] = val;
130 }
131 }
132 } else { // constant Y
133 for (int s = 0; s < p.row_numShapes; s++) {
134 Scalar f = zero;
135 for (int q = 0; q < p.numQuadSub; q++) {
136 f += Vol[q] * S[INDEX2(s,q,p.row_numShapes)];
137 }
138 for (int k = 0; k < p.numEqu; k++) {
139 EM_F[INDEX2(k,s,p.numEqu)] = -f * Y_p[k];
140 EM_F[INDEX2(k,s+p.row_numShapes,p.numEqu)] = f * Y_p[k];
141 }
142 }
143 }
144 }
145 // add the element matrices onto the matrix and
146 // right hand side
147 for (int q = 0; q < p.row_numShapesTotal; q++) {
148 row_index[q] = p.row_DOF[p.elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
149 }
150 if (add_EM_F) {
151 util::addScatter(p.row_numShapesTotal,
152 &row_index[0], p.numEqu, &EM_F[0], F_p,
153 p.row_DOF_UpperBound);
154 }
155 if (add_EM_S)
156 Assemble_addToSystemMatrix(p.S,
157 p.row_numShapesTotal, &row_index[0],
158 p.numEqu, p.col_numShapesTotal,
159 &row_index[0], p.numComp, &EM_S[0]);
160 } // end of isub
161 } // end color check
162 } // end element loop
163 } // end color loop
164 } // end parallel region
165 }
166
167 // instantiate our two supported versions
168 template void Assemble_PDE_System_C<escript::DataTypes::real_t>(
169 const AssembleParameters& p,
170 const escript::Data& D, const escript::Data& Y);
171 template void Assemble_PDE_System_C<escript::DataTypes::cplx_t>(
172 const AssembleParameters& p,
173 const escript::Data& D, const escript::Data& Y);
174
175 } // namespace finley
176

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26