/[escript]/trunk/finley/src/Assemble_PDEMatrix_System2.c
ViewVC logotype

Contents of /trunk/finley/src/Assemble_PDEMatrix_System2.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 616 - (show annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 2 months ago) by elspeth
File MIME type: text/plain
File size: 7829 byte(s)
Copyright added to more source files.

1 /*
2 ************************************************************
3 * Copyright 2006 by ACcESS MNRF *
4 * *
5 * http://www.access.edu.au *
6 * Primary Business: Queensland, Australia *
7 * Licensed under the Open Software License version 3.0 *
8 * http://www.opensource.org/licenses/osl-3.0.php *
9 * *
10 ************************************************************
11 */
12
13 /**************************************************************/
14
15 /* assembles the system of numEq PDEs into the stiffness matrix S: */
16
17 /* -div(A*grad u)-div(B*u)+C*grad u + D*u */
18
19 /* -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m */
20
21 /* u has numComp components. */
22
23 /* Shape of the coefficients: */
24
25 /* A = numEqu x numDim x numComp x numDim */
26 /* B = numDim x numEqu x numComp */
27 /* C = numEqu x numDim x numComp */
28 /* D = numEqu x numComp */
29
30
31 /**************************************************************/
32
33 /* Author: gross@access.edu.au */
34 /* Version: $Id$ */
35
36 /**************************************************************/
37
38
39 #include "Assemble.h"
40
41 /**************************************************************/
42
43 void Finley_Assemble_PDEMatrix_System2(dim_t NS,dim_t numDim,dim_t numQuad,dim_t numEqu,dim_t numComp,
44 double* S,double* DSDX, double* Vol,
45 int NN, double* EM_S,
46 double* A, bool_t extendedA,
47 double* B, bool_t extendedB,
48 double* C, bool_t extendedC,
49 double* D, bool_t extendedD ) {
50 dim_t s,r,k,m,i,j,q;
51 double rtmp;
52
53
54 /**************************************************************/
55 /* process A: */
56 /**************************************************************/
57 if (NULL!=A) {
58 if (extendedA) {
59 for (s=0;s<NS;s++) {
60 for (r=0;r<NS;r++) {
61 for (k=0;k<numEqu;k++) {
62 for (m=0;m<numComp;m++) {
63 for (i=0;i<numDim;i++) {
64 for (j=0;j<numDim;j++) {
65 for (q=0;q<numQuad;q++) {
66 EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=Vol[q]*DSDX[INDEX3(s,i,q,NS,numDim)]*
67 A[INDEX5(k,i,m,j,q,numEqu,numDim,numComp,numDim)]*DSDX[INDEX3(r,j,q,NS,numDim)];
68 }
69 }
70 }
71 }
72 }
73 }
74 }
75 } else {
76 for (s=0;s<NS;s++) {
77 for (r=0;r<NS;r++) {
78 for (i=0;i<numDim;i++) {
79 for (j=0;j<numDim;j++) {
80 rtmp=0;
81 for (q=0;q<numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,i,q,NS,numDim)]*DSDX[INDEX3(r,j,q,NS,numDim)];
82 for (k=0;k<numEqu;k++) {
83 for (m=0;m<numComp;m++) {
84 EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=rtmp*A[INDEX4(k,i,m,j,numEqu,numDim,numComp)];
85 }
86 }
87 }
88 }
89 }
90 }
91 }
92 }
93 /**************************************************************/
94 /* process B: */
95 /**************************************************************/
96 if (NULL!=B) {
97 if (extendedB) {
98 for (s=0;s<NS;s++) {
99 for (r=0;r<NS;r++) {
100 for (k=0;k<numEqu;k++) {
101 for (m=0;m<numComp;m++) {
102 for (i=0;i<numDim;i++) {
103 for (q=0;q<numQuad;q++) {
104 EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=Vol[q]*DSDX[INDEX3(s,i,q,NS,numDim)]*
105 B[INDEX4(k,i,m,q,numEqu,numDim,numComp)]*S[INDEX2(r,q,NS)];
106 }
107 }
108 }
109 }
110 }
111 }
112 } else {
113 for (s=0;s<NS;s++) {
114 for (r=0;r<NS;r++) {
115 for (i=0;i<numDim;i++) {
116 rtmp=0;
117 for (q=0;q<numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,i,q,NS,numDim)]*S[INDEX2(r,q,NS)];
118 for (k=0;k<numEqu;k++) {
119 for (m=0;m<numComp;m++) {
120 EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=rtmp*B[INDEX3(k,i,m,numEqu,numDim)];
121 }
122 }
123 }
124 }
125 }
126 }
127 }
128 /**************************************************************/
129 /* process C: */
130 /**************************************************************/
131 if (NULL!=C) {
132 if (extendedC) {
133 for (s=0;s<NS;s++) {
134 for (r=0;r<NS;r++) {
135 for (k=0;k<numEqu;k++) {
136 for (m=0;m<numComp;m++) {
137 for (j=0;j<numDim;j++) {
138 for (q=0;q<numQuad;q++) {
139 EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=Vol[q]*S[INDEX2(s,q,NS)]*
140 C[INDEX4(k,m,j,q,numEqu,numComp,numDim)]*DSDX[INDEX3(r,j,q,NS,numDim)];
141 }
142 }
143 }
144 }
145 }
146 }
147 } else {
148 for (s=0;s<NS;s++) {
149 for (r=0;r<NS;r++) {
150 for (j=0;j<numDim;j++) {
151 rtmp=0;
152 for (q=0;q<numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,NS)]*DSDX[INDEX3(r,j,q,NS,numDim)];
153 for (k=0;k<numEqu;k++) {
154 for (m=0;m<numComp;m++) {
155 EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=rtmp*C[INDEX3(k,m,j,numEqu,numComp)];
156 }
157 }
158 }
159 }
160 }
161 }
162 }
163 /************************************************************* */
164 /* process D */
165 /**************************************************************/
166 if (NULL!=D) {
167 if (extendedD) {
168 for (s=0;s<NS;s++) {
169 for (r=0;r<NS;r++) {
170 for (k=0;k<numEqu;k++) {
171 for (m=0;m<numComp;m++) {
172 for (q=0;q<numQuad;q++) {
173 EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=Vol[q]*S[INDEX2(s,q,NS)]*
174 D[INDEX3(k,m,q,numEqu,numComp)]*S[INDEX2(r,q,NS)];
175 }
176 }
177 }
178 }
179 }
180 } else {
181 for (s=0;s<NS;s++) {
182 for (r=0;r<NS;r++) {
183 rtmp=0;
184 for (q=0;q<numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,NS)]*S[INDEX2(r,q,NS)];
185 for (k=0;k<numEqu;k++) {
186 for (m=0;m<numComp;m++) {
187 EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=rtmp*D[INDEX2(k,m,numEqu)];
188 }
189 }
190 }
191 }
192 }
193 }
194 }
195
196 /*
197 * $Log$
198 * Revision 1.3 2005/09/15 03:44:21 jgs
199 * Merge of development branch dev-02 back to main trunk on 2005-09-15
200 *
201 * Revision 1.2.2.1 2005/09/07 06:26:17 gross
202 * the solver from finley are put into the standalone package paso now
203 *
204 * Revision 1.2 2005/07/08 04:07:46 jgs
205 * Merge of development branch back to main trunk on 2005-07-08
206 *
207 * Revision 1.1.1.1.2.1 2005/06/29 02:34:47 gross
208 * some changes towards 64 integers in finley
209 *
210 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
211 * initial import of project esys2
212 *
213 * Revision 1.2 2004/07/30 04:37:06 gross
214 * escript and finley are linking now and RecMeshTest.py has been passed
215 *
216 * Revision 1.1 2004/07/02 04:21:13 gross
217 * Finley C code has been included
218 *
219 *
220 */

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26