calcBKBQuad_Lin
| main | Tutorials | Functions | website |
Computes the product B' x K x B for 2D linear quadrilateral elements in vectorized manner This is used internally from calcBKB.
Version : 1.0
Author : George Kourakos
email: giorgk@gmail.com
web : http://groundwater.ucdavis.edu/msim
Date 28-Mar-2014
Department of Land Air and Water
University of California Davis
Contents
Usage
BKB = calcBKBQuad_Lin(B, K, ii)
Input
B: [Nel x N_sh^2] contains the contributions of each element to the final conductance matrix
K : [Nel x Nanis] Hydraulic conductiviy element values. The number of columns is defined by the anisotropy. Maximum number is 3.
ii : In case of nested assembly this indicates the iteration. In vectorized assembly this is not used
Output
BKB: the product B'*K*B
How to compute
In mSim we avoid by hand computations at all costs, therefore we used the symbolic toolbox to perform the vectorized computations. The following code show how we computed the products.
syms b1 b2 b3 b4 syms b5 b6 b7 b8 syms kx ky B=[b1 b2 b3 b4; b5 b6 b7 b8]; BT = [b1 b5; b2 b6; b3 b7; b4 b8]; BKB = BT*[kx 0; 0 ky]*B
BKB = [ kx*b1^2 + ky*b5^2, b1*b2*kx + b5*b6*ky, b1*b3*kx + b5*b7*ky, b1*b4*kx + b5*b8*ky] [ b1*b2*kx + b5*b6*ky, kx*b2^2 + ky*b6^2, b2*b3*kx + b6*b7*ky, b2*b4*kx + b6*b8*ky] [ b1*b3*kx + b5*b7*ky, b2*b3*kx + b6*b7*ky, kx*b3^2 + ky*b7^2, b3*b4*kx + b7*b8*ky] [ b1*b4*kx + b5*b8*ky, b2*b4*kx + b6*b8*ky, b3*b4*kx + b7*b8*ky, kx*b4^2 + ky*b8^2]
Next to convert the above complex expressions to vectorized expresions we used the following loop and copy the output to an editor and with minimum editing we wrote the final expressions without actually write by hand indices and symbols
cnt=0; for i=1:size(B,2) for j=1:size(B,2) cnt=cnt+1; fprintf(['BKB(:,%d) = ' char(BKB(i,j)) ';\n'],cnt); end end
BKB(:,1) = b1^2*kx + b5^2*ky; BKB(:,2) = b1*b2*kx + b5*b6*ky; BKB(:,3) = b1*b3*kx + b5*b7*ky; BKB(:,4) = b1*b4*kx + b5*b8*ky; BKB(:,5) = b1*b2*kx + b5*b6*ky; BKB(:,6) = b2^2*kx + b6^2*ky; BKB(:,7) = b2*b3*kx + b6*b7*ky; BKB(:,8) = b2*b4*kx + b6*b8*ky; BKB(:,9) = b1*b3*kx + b5*b7*ky; BKB(:,10) = b2*b3*kx + b6*b7*ky; BKB(:,11) = b3^2*kx + b7^2*ky; BKB(:,12) = b3*b4*kx + b7*b8*ky; BKB(:,13) = b1*b4*kx + b5*b8*ky; BKB(:,14) = b2*b4*kx + b6*b8*ky; BKB(:,15) = b3*b4*kx + b7*b8*ky; BKB(:,16) = b4^2*kx + b8^2*ky;