/[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 123 - (show annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years, 7 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Assemble_PDEMatrix_System2.c
File MIME type: text/plain
File size: 7105 byte(s)
Merge of development branch back to main trunk on 2005-07-08

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26