Mercurial > repos > public > sbplib

changeset 729:aa8cf3851de8 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update multiblock.DiffOp to work for systems.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 20 Apr 2018 16:56:49 -0700
parents 0aff87f6fb2c
children 9f28cf266f86
files +multiblock/DiffOp.m +scheme/Elastic2dVariable.m
diffstat 2 files changed, 20 insertions(+), 5 deletions(-) [+]
line wrap: on
line diff
diff -r 0aff87f6fb2c -r aa8cf3851de8 +multiblock/DiffOp.m
--- a/+multiblock/DiffOp.m	Fri Apr 20 12:06:38 2018 -0700
+++ b/+multiblock/DiffOp.m	Fri Apr 20 16:56:49 2018 -0700
@@ -53,7 +53,11 @@
 
 
             % Build the differentiation matrix
-            obj.blockmatrixDiv = {grid.Ns, grid.Ns};
+            Ns = zeros(nBlocks,1);
+            for i = 1:nBlocks
+                Ns(i) = length(obj.diffOps{i}.D);
+            end
+            obj.blockmatrixDiv = {Ns, Ns};
             D = blockmatrix.zero(obj.blockmatrixDiv);
             for i = 1:nBlocks
                 D{i,i} = obj.diffOps{i}.D;
diff -r 0aff87f6fb2c -r aa8cf3851de8 +scheme/Elastic2dVariable.m
--- a/+scheme/Elastic2dVariable.m	Fri Apr 20 12:06:38 2018 -0700
+++ b/+scheme/Elastic2dVariable.m	Fri Apr 20 16:56:49 2018 -0700
@@ -64,6 +64,13 @@
 
             h = g.scaling();
             lim = g.lim;
+            if isempty(lim)
+                x = g.x;
+                lim = cell(length(x),1);
+                for i = 1:length(x)
+                    lim{i} = {min(x{i}), max(x{i})};
+                end
+            end
 
             % 1D operators
             ops = cell(dim,1);
@@ -270,6 +277,10 @@
             default_arg('type',{'free','free'});
             default_arg('parameter', []);
 
+            if ~iscell(type)
+                type = {type, type};
+            end
+
             % j is the coordinate direction of the boundary
             j = obj.get_boundary_number(boundary);
             [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
@@ -407,11 +418,11 @@
 
             % Loop over components that penalties end up on
             for i = 1:dim
-                closure = closure - sigma(i,j)*E{i}*RHOi*Hi*e*H_gamma*e'*E{i}';
-                penalty = penalty + sigma(i,j)*E{i}*RHOi*Hi*e*H_gamma*e_v'*E_v{i}';
+                closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}';
+                penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}';
 
-                closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i}';
-                penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i}';
+                closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
+                penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
 
                 % Loop over components that we have interface conditions on
                 for k = 1:dim