comparison +multiblock/DiffOp.m @ 958:72cd29107a9a feature/poroelastic

Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
author Martin Almquist <malmquist@stanford.edu>
date Wed, 05 Dec 2018 18:58:10 -0800
parents 7d4f57725192
children 262b52c3f268
comparison
equal deleted inserted replaced
919:e30aaa4a3e09 958:72cd29107a9a
168 otherwise 168 otherwise
169 error('Unknown boundary indentifier') 169 error('Unknown boundary indentifier')
170 end 170 end
171 end 171 end
172 172
173 % Get a boundary cell of operators, specified by opName for the given boundary/BoundaryGroup
174 function opCell = getBoundaryCellOperator(obj, opName, boundary, blockmatrixDiv)
175 default_arg('blockmatrixDiv', obj.blockmatrixDiv);
176
177 % Get size of cell
178 switch class(boundary)
179 case 'cell'
180 blockId = boundary{1};
181 localCell = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{2});
182 case 'multiblock.BoundaryGroup'
183 blockId = boundary{1}{1};
184 localCell = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{1}{2});
185 otherwise
186 error('Unknown boundary indentifier')
187 end
188
189 % Loop over cell elements and build the boundary operator in each cell
190 opCell = cell(size(localCell));
191 for i = 1:numel(opCell)
192 switch class(boundary)
193 case 'cell'
194 blockId = boundary{1};
195 localOpCell = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{2});
196 localOp = localOpCell{i};
197
198 div = {blockmatrixDiv, size(localOp,2)}
199 blockOp = blockmatrix.zero(div);
200 blockOp{blockId,1} = localOp;
201 op = blockmatrix.toMatrix(blockOp);
202 opCell{i} = op;
203
204 case 'multiblock.BoundaryGroup'
205 op = sparse(size(obj.D,1),0);
206 for j = 1:length(boundary)
207 localCell = obj.getBoundaryCellOperator(opName, boundary{j}, blockmatrixDiv);
208 op = [op, localCell{i}];
209 end
210 opCell{i} = op;
211 otherwise
212 error('Unknown boundary indentifier')
213 end
214 end
215 end
216
173 function op = getBoundaryQuadrature(obj, boundary) 217 function op = getBoundaryQuadrature(obj, boundary)
174 opName = 'H'; 218 opName = 'H';
175 switch class(boundary) 219 switch class(boundary)
176 case 'cell' 220 case 'cell'
177 localOpName = [opName '_' boundary{2}]; 221 localOpName = [opName '_' boundary{2}];