Mercurial > repos > public > sbplib

comparison +scheme/Beam.m @ 220:5df8d20281fe feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Made scheme boundary_condition return a cell array of penalties if there are several of them.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 28 Jun 2016 13:11:14 +0200
parents d095b5396103
children 96dedf892910
comparison
equal deleted inserted replaced
219:f66513508c75 220:5df8d20281fe
60 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 60 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
61 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 61 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
62 % type is a string specifying the type of boundary condition if there are several. 62 % type is a string specifying the type of boundary condition if there are several.
63 % neighbour_scheme is an instance of Scheme that should be interfaced to. 63 % neighbour_scheme is an instance of Scheme that should be interfaced to.
64 % neighbour_boundary is a string specifying which boundary to interface to. 64 % neighbour_boundary is a string specifying which boundary to interface to.
65 function [closure, penalty_e, penalty_d] = boundary_condition(obj,boundary,type) 65 function [closure, penalty] = boundary_condition(obj,boundary,type)
66 default_arg('type','dn'); 66 default_arg('type','dn');
67 67
68 [e, d1, d2, d3, s] = obj.get_boundary_ops(boundary); 68 [e, d1, d2, d3, s] = obj.get_boundary_ops(boundary);
69 gamm = obj.gamm; 69 gamm = obj.gamm;
70 delt = obj.delt; 70 delt = obj.delt;
72 switch type 72 switch type
73 case {'dn'} % Dirichlet-neumann boundary condition 73 case {'dn'} % Dirichlet-neumann boundary condition
74 alpha = obj.alpha; 74 alpha = obj.alpha;
75 75
76 % tau1 < -alpha^2/gamma 76 % tau1 < -alpha^2/gamma
77 % tuning = 2;
77 tuning = 1.1; 78 tuning = 1.1;
78 79
79 tau1 = tuning * alpha/delt; 80 tau1 = tuning * alpha/delt;
80 tau4 = s*alpha; 81 tau4 = s*alpha;
81 82
85 tau = tau1*e+tau4*d3; 86 tau = tau1*e+tau4*d3;
86 sig = sig2*d1+sig3*d2; 87 sig = sig2*d1+sig3*d2;
87 88
88 closure = obj.Hi*(tau*e' + sig*d1'); 89 closure = obj.Hi*(tau*e' + sig*d1');
89 90
90 penalty_e = obj.Hi*tau; 91 penalty{1} = obj.Hi*tau;
91 penalty_d = obj.Hi*sig; 92 penalty{2} = obj.Hi*sig;
92 otherwise % Unknown, boundary condition 93 otherwise % Unknown, boundary condition
93 error('No such boundary condition: type = %s',type); 94 error('No such boundary condition: type = %s',type);
94 end 95 end
95 end 96 end
96 97
104 delt_u = obj.delt; 105 delt_u = obj.delt;
105 106
106 gamm_v = neighbour_scheme.gamm; 107 gamm_v = neighbour_scheme.gamm;
107 delt_v = neighbour_scheme.delt; 108 delt_v = neighbour_scheme.delt;
108 109
109 tuning = 2; 110 % tuning = 2;
111 tuning = 1.1;
110 112
111 alpha_u = obj.alpha; 113 alpha_u = obj.alpha;
112 alpha_v = neighbour_scheme.alpha; 114 alpha_v = neighbour_scheme.alpha;
113 115
114 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning; 116 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning;