Skip to content

Commit

Permalink
add complex coefficients
Browse files Browse the repository at this point in the history
  • Loading branch information
@AndreasStahel
AndreasStahel committed Feb 18, 2025
1 parent 5520617 commit 5862810
Show file tree
Hide file tree
Showing 20 changed files with 2,886 additions and 115 deletions.
29 changes: 15 additions & 14 deletions inst/BVP1DNL.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,8 @@
##@itemize
##@item @var{a} a constant or vector of values of a(x).
##@item @var{a = @@(x)} a function handle to evaluate f(x).
##@item @var{a = @{@@(x,u), @@(x,u)@}} assumes that the function a(x,u) depends on x and u. The two function handles evaluate a(x,u) and the partial derivative a_u(x,u).
##@item @var{a = @{@@(x,u,u'), @@(x,u,u') , @@(x,u,u')@}} assumes that a(x,u,u') depends on x, u and u'. The three function handles evaluate a(x,u,u') and the partial derivatives a_u(x,u,u') and a_u'(x,u,u').
##@item @var{a = @@(x,u)} assumes that the function a(x,u) depends on x and u.
##@item @var{a = @@(x,u,u')} assumes that a(x,u,u') depends on x, u and u'.
##@end itemize
##@item @var{b} constant, vector or function handle to evaluate b(x) at Gauss points
##@item @var{c} constant, vector or function handle to evaluate c(x) at Gauss points
Expand All @@ -53,10 +53,10 @@
##@item @var{BCleft} and @var{BCright} the two boundary conditions
##@itemize
##@item for a Dirichlet condition specify a single value @var{g_D}
##@item for a Neumann condition specify the values @var{[g_N1,g_N2]}
##@item for a Neumann condition specify the two values @var{[g_N1,g_N2]}
##@end itemize
##@item @var{u0} constant, vector or function handle to evaluate u0(x) at the nodes. This is the starting value for the iteration.
##@item @var{options} additional options, given as pairs name/value.
##@item @var{options} additional options, given as pairs name/value.
##@itemize
##@item @var{"tol"} the tolerance for the iteration to stop, given as pair @var{[tolrel,tolabs]} for the relative and absolute tolerance. The iteration stops if the absolute or relative error is smaller than the specified tolerance. RMS (root means square) values are used. If only @var{tolrel} is specified @var{TolAbs=TolRel} is used. The default values are @var{tolrel = tolabs = 1e-5}.
##@item @var{"MaxIter"} the maximal number of iterations to be used. The default value is 10.
Expand Down Expand Up @@ -110,11 +110,11 @@
n = length(x); %% number of nodes

function fun_values = convert2values(fun) %% evaluate at Gauss points
if (length(fun)>1)&&isnumeric(fun) %% a given as vector
if (length(fun)>1)&&isnumeric(fun) %% a given as vector
fun_values = fun;
elseif isnumeric(fun) %% a given as scalar
elseif isnumeric(fun) %% a given as scalar
fun_values = fun*(ones(size(xGauss)));
else %% a given as function handle
else %% a given as function handle
fun_values = fun(xGauss);
endif
endfunction
Expand Down Expand Up @@ -172,6 +172,7 @@
endif %% length(f)==1
endif


%% determine the boundary conditions
if length(BCleft)*length(BCright)==1 %% DD: Dirichlet at both ends
BC = "DD";
Expand Down Expand Up @@ -284,7 +285,7 @@

u_valuesGauss = Nodes2GaussU * u_values;
du_valuesGauss = Nodes2GaussDU * u_values;

if a_NL ; %% update the coefficient a, if necessary
switch nargin(a)
case 2 %% a depends on x and u
Expand All @@ -293,7 +294,7 @@
a_values = a(xGauss,u_valuesGauss,du_valuesGauss);
endswitch
endif % a_NL

if f_NL %%% evaluate f and apply one Newton step
switch length(f) %% evaluate f and its partial derivatives
case 1 %% f(x)
Expand All @@ -304,7 +305,7 @@
dfdu_valuesGauss = f{2}(xGauss,u_valuesGauss); %% df/du(x,u) at Gauss
A2 = GenerateFEM1D(interval,a_values,b,c_values-d_values.*dfdu_valuesGauss,d);
case 3 %% f(x,u,u')
du_values = FEM1DEvaluateDu(x,u_values); %% u' at nodes
du_values = FEM1DEvaluateDu(x,u_values); %% u' at nodes
f_values = f{1}(x,u_values,du_values); %% f(x,u,u') at nodes
dfdu_valuesGauss = f{2}(xGauss,u_valuesGauss,du_valuesGauss);% df/du
dfdupr_valuesGauss = f{3}(xGauss,u_valuesGauss,du_valuesGauss);% df/du'
Expand All @@ -319,7 +320,7 @@
else
Au = A*u_values;
endif

switch BC %% solve the problem for phi and make one Newton step
case "DD"
phi = SolveSystemPhi(A2,-Au(2:end-1)+M*f_values,BCleft,BCright);
Expand All @@ -336,7 +337,7 @@
case 2 %% f(x,u)
f_values = f{1}(x,u_values); %% f(x,u) at nodes
case 3 %% f(x,u,u')
du_values = FEM1DEvaluateDu(x,u_values); %% u' at nodes
du_values = FEM1DEvaluateDu(x,u_values); %% u' at nodes
f_values = f{1}(x,u_values,du_values); %% f(x,u,u') at nodes
endswitch %% length(f)
endif %% f_NL
Expand All @@ -355,7 +356,7 @@
case 2 %% f(x,u)
f_values = f{1}(x,u_values); %% f(x,u) at nodes
case 3 %% f(x,u,u')
du_values = FEM1DEvaluateDu(x,u_values); %% u' at nodes
du_values = FEM1DEvaluateDu(x,u_values); %% u' at nodes
f_values = f{1}(x,u_values,du_values); %% f(x,u,u') at nodes
endswitch %% length(f)
endif %% f_NL
Expand All @@ -379,7 +380,7 @@
endif %% strcmp

u = u_values; %% assign the return value u

if f_NL
AbsError2 = mean([AbsError,norm(phi)]);
else
Expand Down
94 changes: 77 additions & 17 deletions inst/BVP2D.m
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
## Copyright (C) 2020 Andreas Stahel
## Copyright (C) 2025 Andreas Stahel
##
## This program is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
Expand All @@ -16,11 +16,11 @@


## Author: Andreas Stahel <andreas.stahel@gmx.com>
## Created: 2020-03-30
## Created: 2025-01-04


## -*- texinfo -*-
## @deftypefn{function file}{}@var{u} = BVP2D(@var{mesh},@var{a},@var{b0},@var{bx},@var{by},@var{f},@var{gD},@var{gN1},@var{gN2})
## @deftypefn{function file }{}@var{ u} = BVP2D(@var{mesh},@var{a},@var{b0},@var{bx},@var{by},@var{f},@var{gD},@var{gN1},@var{gN2},@var{options})
##
## Solve an elliptic boundary value problem
##
Expand All @@ -42,7 +42,13 @@
##A=[axx,axy;axy,ayy] given by the row vector [axx,ayy,axy].
##It can be given as row vector or as string with the function name or
##as nx3 matrix with the values at the Gauss points.

##@item @var{options} additional options, given as pairs name/value.
##Currently only real or complex coefficient problems can be selected by
##@var{"TYPE"} and the possible values are
##@itemize
##@item @var{"real"}: all coefficients are real (default)
##@item @var{"complex"}: some coefficients might be complex
##@end itemize
##@end itemize
##
##return value
Expand All @@ -58,20 +64,74 @@
## @c END_CUT_TEXINFO
## @end deftypefn

function u = BVP2D(Mesh,a,b0,bx,by,f,gD,gN1,gN2)
if nargin ~=9
print_usage();
endif
function u = BVP2D(Mesh,a,b0,bx,by,f,gD,gN1,gN2,varargin)
if nargin < 9
print_usage();
endif

Type = 'REAL'; %% default value
if (~isempty(varargin))
for cc = 1:2:length(varargin)
switch toupper(varargin{cc})
case {'TYPE'}
Type = toupper(varargin{cc+1});
otherwise
error('Invalid optional argument, %s. Possible value TYPE',varargin{cc});
endswitch % switch
endfor % for
endif % if isempty

switch Mesh.type
case 'linear' %% first order elements
[A,b] = FEMEquation (Mesh,a,b0,bx,by,f,gD,gN1,gN2);
case 'quadratic' %% second order elements
[A,b] = FEMEquationQuad(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
case 'cubic' %% third order elements
[A,b] = FEMEquationCubic(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
endswitch
if ((strcmp(Type,'REAL')==0)&&(strcmp(Type,'COMPLEX')==0))
error('wrong TYPE, possible values are REAL or COMPLEX')
endif

switch Mesh.type
case 'linear' %% first order elements
switch Type
case 'REAL'
if exist("FEMEquation.oct")==3
[A,b] = FEMEquation(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
else
[A,b] = FEMEquationM(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
endif
case 'COMPLEX'
if exist("FEMEquationComplex.oct")==3
[A,b] = FEMEquationComplex(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
else
[A,b] = FEMEquationM(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
endif
endswitch
case 'quadratic' %% second order elements
switch Type
case 'REAL'
if exist("FEMEquationQuad.oct")==3
[A,b] = FEMEquationQuad(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
else
[A,b] = FEMEquationQuadM(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
endif
case 'COMPLEX'
if exist("FEMEquationQuadComplex.oct")==3
[A,b] = FEMEquationQuadComplex(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
else
[A,b] = FEMEquationQuadM(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
endif
endswitch
case 'cubic' %% third order elements
switch Type
case 'REAL'
if exist("FEMEquationCubic.oct")==3
[A,b] = FEMEquationCubic(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
else
[A,b] = FEMEquationCubicM(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
endif
case 'COMPLEX'
if exist("FEMEquationCubicComplex.oct")==3
[A,b] = FEMEquationCubicComplex(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
else
[A,b] = FEMEquationCubicM(Mesh,a,b0,bx,by,f,gD,gN1,gN2);
endif
endswitch
endswitch

u = FEMSolve(Mesh,A,b,gD); %% solve the resulting system
u = FEMSolve(Mesh,A,b,gD); %% solve the resulting system
endfunction
Loading

0 comments on commit 5862810

Please sign in to comment.