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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4442 - (show annotations)
Fri Jun 7 07:04:44 2013 UTC (6 years, 1 month ago) by caltinay
File size: 8233 byte(s)
cosmetic changes.

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

Properties

Name Value
svn:mergeinfo /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 /release/3.0/finley/src/Assemble_PDE_System2_C.cpp:2591-2601 /trunk/finley/src/Assemble_PDE_System2_C.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26