Skip to content

Commit

Permalink
add anisotropic code
Browse files Browse the repository at this point in the history
  • Loading branch information
@AndreasStahel
AndreasStahel committed Nov 19, 2024
1 parent 3485713 commit aea35e5
Show file tree
Hide file tree
Showing 9 changed files with 705 additions and 463 deletions.
5 changes: 5 additions & 0 deletions inst/BVP2D.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,11 @@
##@*Any constant function can be given by its scalar value.
##@*The functions @var{a},@var{b0},@var{bx},@var{by} and @var{f} may also be given as vectors
##with the values of the function at the Gauss points.
##@item The coefficient @var{a} can also be a symmetric matrix
##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.

##@end itemize
##
##return value
Expand Down
4 changes: 4 additions & 0 deletions inst/BVP2Deig.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,10 @@
##@*Any constant function can be given by its scalar value.
##@*The functions @var{a},@var{b0} and @var{w} may also be given as vectors
##with the values of the function at the Gauss points.
##@item The coefficient @var{a} can also be a symmetric matrix
## 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{nVec} number of smallest eigenvalues to be computed
##@item@var{tol} optional tolerance for the eigenvalue iteration, default 1e-5
##@end itemize
Expand Down
4 changes: 4 additions & 0 deletions inst/BVP2Dsym.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,10 @@
##@*Any constant function can be given by its scalar value.
##@*The functions @var{a},@var{b0} and @var{f} may also be given as vectors
##with the values of the function at the Gauss points.
##@item The coefficient @var{a} can also be a symmetric matrix
## 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.
##@end itemize
##
##return value
Expand Down
48 changes: 42 additions & 6 deletions inst/FEMEquation.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,11 @@
% 'a','b0','bx','by','f','gD','gN1','gN2' are the functions and coefficients
% for the boundary value problem. They can be given as a scalar value
% or as a string with the function name
% The coefficient 'a' can also be a symmetric matrix 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.
%
%
% -div(a*grad u-u*(bx,by)) + b0*u = f in domain
% u = gD on Dirichlet section of the boundary
Expand All @@ -22,12 +27,37 @@

nElem = size(Mesh.elem,1); nGP = size(Mesh.GP,1) ;

if ischar(aFunc)
aV = reshape(feval(aFunc,Mesh.GP,Mesh.GPT),nGP/nElem,nElem);
elseif isscalar(aFunc)
isotropic = 1; %% flag to mark that the coefficient a is a scalar

if ischar(aFunc) %% given as string
aV = feval(aFunc,Mesh.GP,Mesh.GPT);
sizeA = size(aV);
if sizeA(2)==3
isotropic = 0;
a11V = reshape(aV(:,1),nGP/nElem,nElem);
a22V = reshape(aV(:,2),nGP/nElem,nElem);
a12V = reshape(aV(:,3),nGP/nElem,nElem);
else
aV = reshape(aV,nGP/nElem,nElem);
endif
elseif isscalar(aFunc) %% given as one scalar value
aV = aFunc*ones(nGP/nElem,nElem);
else
aV = reshape(aFunc,nGP/nElem,nElem);
else %% given as values at the Gauss points
aSize = size(aFunc);
if aSize(2)==3 %% three coefficients, values at Gauss points
isotropic = 0;
if (aSize(1)==1) %% three scalar values
a11V = aFunc(1)*ones(nGP/nElem,nElem);
a22V = aFunc(2)*ones(nGP/nElem,nElem);
a12V = aFunc(3)*ones(nGP/nElem,nElem);
else %% matrix with values at all Gauss points
a11V = reshape(aFunc(:,1),nGP/nElem,nElem);
a22V = reshape(aFunc(:,2),nGP/nElem,nElem);
a12V = reshape(aFunc(:,3),nGP/nElem,nElem);
endif
else %% one coefficient, values at Gauss points
aV = reshape(aFunc,nGP/nElem,nElem);
endif
endif

if ischar(b0Func)
Expand Down Expand Up @@ -77,7 +107,13 @@
G = [cor(3,2)-cor(2,2),cor(1,2)-cor(3,2),cor(2,2)-cor(1,2);...
cor(2,1)-cor(3,1),cor(3,1)-cor(1,1),cor(1,1)-cor(2,1)];
Mdiffusion = (bxV(:,k)*G(1,:) + byV(:,k)*G(2,:))'*M/6;
mat = sum(aV(:,k))/(12*area)*G'*G + area/3*M*diag(bV(:,k))*M+Mdiffusion;
if isotropic
mat = sum(aV(:,k))/(12*area)*G'*G + area/3*M*diag(bV(:,k))*M+Mdiffusion;
else
a11 = sum(a11V(:,k));a22 = sum(a22V(:,k));a12 = sum(a12V(:,k));
a = [a11,a12;a12,a22];
mat = 1/(12*area)*G'*a*G + area/3*M*diag(bV(:,k))*M+Mdiffusion;
endif
vec = -area/3*M*fV(:,k);
dofs = Mesh.node2DOF(Mesh.elem(k,:));
for k1 = 1:3
Expand Down
57 changes: 47 additions & 10 deletions inst/FEMEquationCubic.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,10 @@
% 'a','b0','bx',by','f','gD','gN1','gN2' are the functions and coefficients
% for the boundary value problem. They can be given as a scalar value
% or as a string with the function name
% The coefficient 'a' can also be a symmetric matrix 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.
%
% -div(a*grad u-u*(bx,by)) + b0*u = f in domain%
% u = gD on Dirichlet section of the boundary
Expand All @@ -22,12 +26,37 @@
nElem = size(Mesh.elem,1);
nGP = size(Mesh.GP,1) ;

if ischar(aFunc)
aV = reshape(feval(aFunc,Mesh.GP,Mesh.nodesT),nGP/nElem,nElem);
elseif isscalar(aFunc)
isotropic = 1; %% flag to mark that the coefficient a is a scalar

if ischar(aFunc) %% given as string
aV = feval(aFunc,Mesh.GP,Mesh.GPT);
sizeA = size(aV);
if sizeA(2)==3
isotropic = 0;
a11V = reshape(aV(:,1),nGP/nElem,nElem);
a22V = reshape(aV(:,2),nGP/nElem,nElem);
a12V = reshape(aV(:,3),nGP/nElem,nElem);
else
aV = reshape(aV,nGP/nElem,nElem);
endif
elseif isscalar(aFunc) %% given as one scalar value
aV = aFunc*ones(nGP/nElem,nElem);
else
aV = reshape(aFunc,nGP/nElem,nElem);
else %% given as values at the Gauss points
aSize = size(aFunc);
if aSize(2)==3 %% three coefficients, values at Gauss points
isotropic = 0;
if (aSize(1)==1) %% three scalar values
a11V = aFunc(1)*ones(nGP/nElem,nElem);
a22V = aFunc(2)*ones(nGP/nElem,nElem);
a12V = aFunc(3)*ones(nGP/nElem,nElem);
else %% matrix with values at all Gauss points
a11V = reshape(aFunc(:,1),nGP/nElem,nElem);
a22V = reshape(aFunc(:,2),nGP/nElem,nElem);
a12V = reshape(aFunc(:,3),nGP/nElem,nElem);
endif
else %% one coefficient, values at Gauss points
aV = reshape(aFunc,nGP/nElem,nElem);
endif
endif

if ischar(b0Func)
Expand Down Expand Up @@ -117,13 +146,21 @@
G = [cor(3,2)-cor(2,2),cor(1,2)-cor(3,2),cor(2,2)-cor(1,2);...
cor(2,1)-cor(3,1),cor(3,1)-cor(1,1),cor(1,1)-cor(2,1)];
B = M'*diag(w.*bV(:,k))*M;
Mtemp = -G(1,2)*Mxi-G(1,3)*Mnu;
Ax = Mtemp'*diag(w.*aV(:,k))*Mtemp;
Mtemp = -G(2,2)*Mxi-G(2,3)*Mnu;
Ay = Mtemp'*diag(w.*aV(:,k))*Mtemp;
MtempX = -G(1,2)*Mxi-G(1,3)*Mnu;
MtempY = -G(2,2)*Mxi-G(2,3)*Mnu;
if isotropic
A11 = MtempX'*diag(w.*aV(:,k))*MtempX;
A22 = MtempY'*diag(w.*aV(:,k))*MtempY;
mat = 0.5/area*(A11+A22) + 2*area*B;
else
A11 = MtempX'*diag(w.*a11V(:,k))*MtempX;
A22 = MtempY'*diag(w.*a22V(:,k))*MtempY;
A12 = MtempY'*diag(w.*a12V(:,k))*MtempX;
mat = 0.5/area*(A11+A22+2*A12) + 2*area*B;
endif
Abxy = ((G(1,2)*Mxi + G(1,3)*Mnu)'*diag(w.*bxV(:,k))...
+(G(2,2)*Mxi + G(2,3)*Mnu)'*diag(w.*byV(:,k)) )*M;
mat = 0.5/area*(Ax+Ay) + 2*area*B + Abxy;
mat = mat + Abxy;
vec = -2*area*M'*(w.*fV(:,k));
dofs = Mesh.node2DOF(Mesh.elem(k,:));
for k1 = 1:10
Expand Down
57 changes: 47 additions & 10 deletions inst/FEMEquationQuad.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,10 @@
% 'a','b0', 'bx','by','f','gD','gN1','gN2' are the functions and coefficients
% for the boundary value problem. They can be given as a scalar value
% or as a string with the function name
% The coefficient 'a' can also be a symmetric matrix 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.
%
% -div(a*grad u-u*(bx,by)) + b0*u = f in domain%
% u = gD on Dirichlet section of the boundary
Expand All @@ -23,12 +27,37 @@
nElem = size(Mesh.elem,1);
nGP = size(Mesh.GP,1) ;

if ischar(aFunc)
aV = reshape(feval(aFunc,Mesh.GP,Mesh.GPT),nGP/nElem,nElem);
elseif isscalar(aFunc)
isotropic = 1; %% flag to mark that the coefficient a is a scalar

if ischar(aFunc) %% given as string
aV = feval(aFunc,Mesh.GP,Mesh.GPT);
sizeA = size(aV);
if sizeA(2)==3
isotropic = 0;
a11V = reshape(aV(:,1),nGP/nElem,nElem);
a22V = reshape(aV(:,2),nGP/nElem,nElem);
a12V = reshape(aV(:,3),nGP/nElem,nElem);
else
aV = reshape(aV,nGP/nElem,nElem);
endif
elseif isscalar(aFunc) %% given as one scalar value
aV = aFunc*ones(nGP/nElem,nElem);
else
aV = reshape(aFunc,nGP/nElem,nElem);
else %% given as values at the Gauss points
aSize = size(aFunc);
if aSize(2)==3 %% three coefficients, values at Gauss points
isotropic = 0;
if (aSize(1)==1) %% three scalar values
a11V = aFunc(1)*ones(nGP/nElem,nElem);
a22V = aFunc(2)*ones(nGP/nElem,nElem);
a12V = aFunc(3)*ones(nGP/nElem,nElem);
else %% matrix with values at all Gauss points
a11V = reshape(aFunc(:,1),nGP/nElem,nElem);
a22V = reshape(aFunc(:,2),nGP/nElem,nElem);
a12V = reshape(aFunc(:,3),nGP/nElem,nElem);
endif
else %% one coefficient, values at Gauss points
aV = reshape(aFunc,nGP/nElem,nElem);
endif
endif

if ischar(b0Func)
Expand Down Expand Up @@ -88,13 +117,21 @@
G = [cor(3,2)-cor(2,2),cor(1,2)-cor(3,2),cor(2,2)-cor(1,2);...
cor(2,1)-cor(3,1),cor(3,1)-cor(1,1),cor(1,1)-cor(2,1)];
B = M'*diag(w.*bV(:,k))*M;
Mtemp = -G(1,2)*Mxi-G(1,3)*Mnu;
Ax = Mtemp'*diag(w.*aV(:,k))*Mtemp;
Mtemp = -G(2,2)*Mxi-G(2,3)*Mnu;
Ay = Mtemp'*diag(w.*aV(:,k))*Mtemp;
MtempX = -G(1,2)*Mxi-G(1,3)*Mnu;
MtempY = -G(2,2)*Mxi-G(2,3)*Mnu;
if isotropic
A11 = MtempX'*diag(w.*aV(:,k))*MtempX;
A22 = MtempY'*diag(w.*aV(:,k))*MtempY;
mat = 0.5/area*(A11+A22) + 2*area*B;
else
A11 = MtempX'*diag(w.*a11V(:,k))*MtempX;
A22 = MtempY'*diag(w.*a22V(:,k))*MtempY;
A12 = MtempY'*diag(w.*a12V(:,k))*MtempX;
mat = 0.5/area*(A11+A22+2*A12) + 2*area*B;
endif
Abxy = ((G(1,2)*Mxi + G(1,3)*Mnu)'*diag(w.*bxV(:,k))...
+(G(2,2)*Mxi + G(2,3)*Mnu)'*diag(w.*byV(:,k)) )*M;
mat = 0.5/area*(Ax+Ay) + 2*area*B + Abxy;
mat = mat + Abxy;
vec = -2*area*M'*(w.*fV(:,k));
dofs = Mesh.node2DOF(Mesh.elem(k,:));
for k1 = 1:6
Expand Down
Loading

0 comments on commit aea35e5

Please sign in to comment.