Skip to content

Commit

Permalink
Merge remote-tracking branch 'private-repo/Cahn-Hilliard' into merge_…
Browse files Browse the repository at this point in the history
…from_private
  • Loading branch information
rafavzqz committed Jul 29, 2024
2 parents 6e16996 + c94fe7c commit a233bc4
Show file tree
Hide file tree
Showing 137 changed files with 10,397 additions and 69 deletions.
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -9,3 +9,4 @@
*.log
*.out
*.swp
*.asv
2 changes: 1 addition & 1 deletion geopdes/CITATION
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ The BibTeX entries for LaTeX users are:
@article{geopdes,
Author = {C. de Falco and A. Reali and R. V{\'a}zquez},
Journal = {Advances in Engineering Software},
Number = {12},
Number = {3},
Pages = {1020-1034},
Title = {Geo{PDE}s: A research tool for Isogeometric Analysis of {PDE}s},
Volume = {42},
Expand Down
6 changes: 3 additions & 3 deletions geopdes/DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
Name: geopdes
Version: 3.2.2
Date: 2021-04-07
Version: 3.4.0
Date: 2024-07-24
Author: Rafael Vazquez, Carlo de Falco
Maintainer: Rafael Vazquez
Title: Iso-Geometric Analysis
Description: A research tool for Iso-Geometric Analysis of PDEs
Depends: octave (>= 5.1), nurbs (>= 1.4.1)
Depends: octave (>= 6.4), nurbs (>= 1.4.1)
License: GPLv3
36 changes: 35 additions & 1 deletion geopdes/INDEX
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,19 @@ Solve: solvers for some simple problems, that you can modify to solve your own p
solve_bilaplace_gradgrad_2d_iso
solve_adv_diff_2d
solve_linear_elasticity
solve_linear_elasticity_scalarspace
solve_maxwell_eig
solve_maxwell_eig_mixed1
solve_maxwell_eig_mixed2_2d
solve_maxwell_src
solve_navier_stokes
solve_stokes
solve_kirchhoff_love_shell
solve_kirchhoff_love_shell_scalarspace
solve_BE_Beam
solve_cahn_hilliard
cahn_hilliard/generalized_alpha_step_cahn_hilliard
cahn_hilliard/Res_K_cahn_hilliard
Geometry: read geometry files, that can be created with the NURBS package
geo_nurbs
geo_load
Expand Down Expand Up @@ -73,6 +78,8 @@ Space: functions to generate and compute discrete spline and NURBS spaces
@sp_scalar/sp_to_vtk
@sp_scalar/sp_drchlt_l2_proj
@sp_scalar/sp_weak_drchlt_bc_laplace
@sp_scalar/op_penalty_dudn
@sp_scalar/op_nitsche_consistency_cahn_hilliard
@sp_scalar/sp_get_basis_functions
@sp_scalar/sp_get_cells
@sp_scalar/sp_get_neighbors
Expand Down Expand Up @@ -159,6 +166,7 @@ Multipatch: functions of different kinds to solve problems in domains defined by
mp_solve_maxwell_src
mp_solve_stokes
mp_solve_stokes_div_conforming
mp_solve_cahn_hilliard_C1
msh_restrict_to_patches
msh_multipatch
@msh_multipatch/msh_evaluate_element_list
Expand All @@ -180,11 +188,37 @@ Multipatch: functions of different kinds to solve problems in domains defined by
@sp_multipatch/sp_weak_drchlt_bc_laplace
@sp_multipatch/sp_weak_drchlt_bc_stokes
@sp_multipatch/sp_refine
sp_multipatch_C1
@sp_multipatch_C1/sp_l2_error
@sp_multipatch_C1/sp_h1_error
@sp_multipatch_C1/sp_h2_error
@sp_multipatch_C1/sp_h2_equiv_lap_error
@sp_multipatch_C1/sp_plot_solution
@sp_multipatch_C1/sp_to_vtk
@sp_multipatch_C1/sp_eval_phys
@sp_multipatch_C1/sp_evaluate_element_list
@sp_multipatch_C1/sp_get_cells
@sp_multipatch_C1/sp_get_basis_functions
@sp_multipatch_C1/sp_get_functions_on_patch
@sp_multipatch_C1/sp_get_local_interior_functions
@sp_multipatch_C1/sp_get_neighbors
@sp_multipatch_C1/sp_compute_Cpatch
@sp_multipatch_C1/sp_compute_Cpatch_vector
@sp_multipatch_C1/sp_bilaplacian_drchlt_C1_exact
@sp_multipatch_C1/sp_bilaplacian_drchlt_C1
@sp_multipatch_C1/sp_drchlt_C1_shells
@sp_multipatch_C1/sp_nitsche_KL_rotation
@sp_multipatch_C1/sp_weak_drchlt_bc_laplace
@sp_multipatch_C1/op_penalty_dudn
@sp_multipatch_C1/op_nitsche_consistency_cahn_hilliard
@sp_multipatch_C1/sp_refine
Utilities: other functions in the package
collocation_csp
knt_derham
kntunclamp
msh_to_vtk
newtons_method
get_boundary_indices
get_volumetric_indices
get_volumetric_indices
isotropic_stiffness
op_KL_rotation_fluxes
7 changes: 7 additions & 0 deletions geopdes/NEWS
Original file line number Diff line number Diff line change
@@ -1,3 +1,10 @@
Summary of important changes for geopdes-3.4.0:
-----------------------------------------------------------------------------
* Added new class SP_MULTIPATCH_C1.
* Added solvers for C1 spaces: Laplace, bilaplace, linear elast., Kirchhoff-Love shell.
* Added solver for linear elasticity and Kirchhoff-Love using scalar spaces.
* Added solver for Cahn-Hilliard equations, in single patch and multipatch.

Summary of important changes for geopdes-3.3.0:
-----------------------------------------------------------------------------
* Added functions to compute the exterior derivative.
Expand Down
4 changes: 2 additions & 2 deletions geopdes/PKG_ADD
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,8 @@ dirlist = {"examples", "examples/geometry_files", "examples/geometry_file
"examples/elasticity", "examples/elasticity/data_files", ...
"examples/fluid", "examples/fluid/data_files", ...
"examples/maxwell", "examples/maxwell/data_files", ...
"examples/kirchhoff_love_shell", ...
"geometry", "utils", "msh", "obsolete", "operators", "solve", "space", "multipatch"};
"examples/kirchhoff_love_shell", "examples/cahn-hilliard", ...
"geometry", "utils", "msh", "obsolete", "operators", "solve", "space", "multipatch", "solve/cahn_hilliard"};

dir = fileparts (mfilename ("fullpath"));

Expand Down
4 changes: 2 additions & 2 deletions geopdes/PKG_DEL
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,8 @@ dirlist = {"examples", "examples/geometry_files", "examples/geometry_file
"examples/elasticity", "examples/elasticity/data_files", ...
"examples/fluid", "examples/fluid/data_files", ...
"examples/maxwell", "examples/maxwell/data_files", ...
"examples/kirchhoff_love_shell", ...
"geometry", "utils", "msh", "obsolete", "operators", "solve", "space", "multipatch"};
"examples/kirchhoff_love_shell", "examples/cahn-hilliard", ...
"geometry", "utils", "msh", "obsolete", "operators", "solve", "space", "multipatch", "solve/cahn_hilliard"};

[basename,dir] = fileparts(fileparts (mfilename ("fullpath")));

Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
% TEST_PLATE_MIXED_BC_G_NMNN: data function for Neumann boundary condition.

function g = bilaplacian_Lshaped_g_nmnn_8patches (x, y, ind)
[th, r] = cart2pol (x, y);
th = (th < 0).*(2*acos(-1) + th) + (th >= 0) .* th;
% f = @(z) sin(3*pi/2*z).^2 - z.^2 * sin(3*pi/2)^2;
z = 0.544483736782464; %z = fsolve(f, 0.5) The one on the left is more precise, computed with bisection method

C1 = 1/(z-1)*sin(3*pi/2*(z-1)) - 1/(z+1)*sin(3*pi/2*(z+1));
C2 = cos(3*pi/2*(z-1)) - cos(3*pi/2*(z+1));
F1 = cos((z-1)*th) - cos((z+1)*th);
F2 = 1/(z-1)*sin((z-1)*th) - 1/(z+1)*sin((z+1)*th);
DF1 = -(z-1)*sin((z-1)*th)+(z+1)*sin((z+1)*th);
DF2 = cos((z-1)*th)-cos((z+1)*th);

psi = C1*F1 - C2*F2; %C1 and C2 are constants
Dpsi = C1*DF1 - C2*DF2;

switch (ind)
case {4} %n=(-1,0)
g = -((z+1)*r.^(z).* psi.*cos(th) - r.^z.*Dpsi.*sin(th));
case {2,6} %n=(1,0)
g = (z+1)*r.^(z).* psi.*cos(th) - r.^z.*Dpsi.*sin(th);
case {1,3} %n=(0,-1)
g = -((z+1)*r.^(z).* psi.*sin(th) + r.^z.*Dpsi.*cos(th));
case {5} %n=(0,1)
g = (z+1)*r.^(z).* psi.*sin(th) + r.^z.*Dpsi.*cos(th);

otherwise
error ('g_nmnn: unknown reference number')
end

end

21 changes: 21 additions & 0 deletions geopdes/inst/examples/base/data_files/bilaplacian_rhs_Lshaped.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
function uex = bilaplacian_rhs_Lshaped (x, y)
[th, r] = cart2pol (x, y);
th = (th < 0).*(2*acos(-1) + th) + (th >= 0) .* th;
z = 0.544483736782464;

C1 = 1/(z-1)*sin(3*pi/2*(z-1)) - 1/(z+1)*sin(3*pi/2*(z+1));
C2 = cos(3*pi/2*(z-1)) - cos(3*pi/2*(z+1));
F1 = cos((z-1)*th) - cos((z+1)*th);
F1_2 = -(z-1)^2*cos((z-1)*th) + (z+1)^2*cos((z+1)*th);
F1_4 = (z-1)^4*cos((z-1)*th) - (z+1)^4*cos((z+1)*th);
F2 = 1/(z-1)*sin((z-1)*th) - 1/(z+1)*sin((z+1)*th);
F2_2 = -(z-1)*sin((z-1)*th) + (z+1)*sin((z+1)*th);
F2_4 = (z-1)^3*sin((z-1)*th) - (z+1)^3*sin((z+1)*th);

psi = C1*F1 - C2*F2;
psi2 = C1*F1_2 - C2*F2_2;
psi4 = C1*F1_4 - C2*F2_4;

uex = (r.^(z-3)).*(((z-1)^2)*(((z+1)^2)*psi+psi2)+((z+1)^2)*psi2+psi4);

end
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
function uex = solution_bilaplacian_Lshaped (x, y)
[th, r] = cart2pol (x, y);
th = (th < 0).*(2*acos(-1) + th) + (th >= 0) .* th;
% f = @(z) sin(3*pi/2*z).^2 - z.^2 * sin(3*pi/2)^2;
z = 0.544483736782464; %z = fsolve(f, 0.5) The one on the left is more precise, computed with bisection method

C1 = 1/(z-1)*sin(3*pi/2*(z-1)) - 1/(z+1)*sin(3*pi/2*(z+1));
C2 = cos(3*pi/2*(z-1)) - cos(3*pi/2*(z+1));
F1 = cos((z-1)*th) - cos((z+1)*th);
F2 = 1/(z-1)*sin((z-1)*th) - 1/(z+1)*sin((z+1)*th);

psi = C1*F1 - C2*F2; %C1 and C2 are constants

uex = r.^(z+1).* psi;
end

Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
function graduex = solution_bilaplacian_Lshaped_grad (x, y)
graduex = zeros(2, size(x,1), size(x,2));
[th, r] = cart2pol (x, y);
th = (th < 0).*(2*acos(-1) + th) + (th >= 0) .* th;
% f = @(z) sin(3*pi/2*z).^2 - z.^2 * sin(3*pi/2)^2;
z = 0.544483736782464; %z = fsolve(f, 0.5) The one on the left is more precise, computed with bisection method

C1 = 1/(z-1)*sin(3*pi/2*(z-1)) - 1/(z+1)*sin(3*pi/2*(z+1));
C2 = cos(3*pi/2*(z-1)) - cos(3*pi/2*(z+1));
F1 = cos((z-1)*th) - cos((z+1)*th);
F2 = 1/(z-1)*sin((z-1)*th) - 1/(z+1)*sin((z+1)*th);

psi = C1*F1 - C2*F2; %C1 and C2 are constants
Dpsi = C1*(-(z-1)*sin((z-1)*th)+(z+1)*sin((z+1)*th)) - C2*(cos((z-1)*th)-cos((z+1)*th));

graduex(1,:,:) = (z+1)*r.^(z).* psi.*cos(th) - r.^z.*Dpsi.*sin(th);
graduex(2,:,:) = (z+1)*r.^(z).* psi.*sin(th) + r.^z.*Dpsi.*cos(th);
end

Original file line number Diff line number Diff line change
@@ -0,0 +1,36 @@
function hessian = solution_bilaplacian_Lshaped_hessian (x, y)
hessian = zeros(2,2,size(x,1), size(x,2));
[th, r] = cart2pol (x, y);
th = (th < 0).*(2*acos(-1) + th) + (th >= 0) .* th;
z = 0.544483736782464;

C1 = 1/(z-1)*sin(3*pi/2*(z-1)) - 1/(z+1)*sin(3*pi/2*(z+1));
C2 = cos(3*pi/2*(z-1)) - cos(3*pi/2*(z+1));
F1 = cos((z-1)*th) - cos((z+1)*th);
F2 = 1/(z-1)*sin((z-1)*th) - 1/(z+1)*sin((z+1)*th);

psi = C1*F1 - C2*F2; %C1 and C2 are constants
Dpsi = C1*(-(z-1)*sin((z-1)*th)+(z+1)*sin((z+1)*th)) - C2*(cos((z-1)*th)-cos((z+1)*th));
D2psi = C1*(-((z-1)^2)*cos((z-1)*th)+((z+1)^2)*cos((z+1)*th)) - C2*(-(z-1)*sin((z-1)*th)+(z+1)*sin((z+1)*th));

ur=(z+1)*r.^(z).*psi;
uth=r.^(z+1).*Dpsi;
urr=z*(z+1)*r.^(z-1).*psi;
uthr=(z+1)*r.^(z).*Dpsi;
uthth=r.^(z+1).*D2psi;

uxr=cos(th).*urr-sin(th).*(r.*uthr-uth)./(r.^2);
uxth=-sin(th).*ur+cos(th).*uthr-cos(th).*uth./r-sin(th).*uthth./r;
uyr=sin(th).*urr+cos(th).*(r.*uthr-uth)./(r.^2);
uyth=cos(th).*ur+sin(th).*uthr-sin(th).*uth./r+cos(th).*uthth./r;

uxx=cos(th).*uxr-sin(th).*uxth./r;
uxy=sin(th).*uxr+cos(th).*uxth./r;
uyy=sin(th).*uyr+cos(th).*uyth./r;

hessian(1,1,:,:) = uxx;
hessian(1,2,:,:) = uxy;
hessian(2,1,:,:) = uxy;
hessian(2,2,:,:) = uyy;

end
48 changes: 48 additions & 0 deletions geopdes/inst/examples/base/ex_bilaplacian_Lshaped_8patch.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
% PHYSICAL DATA OF THE PROBLEM
clear problem_data
problem_data.geo_name = 'geo_Lshaped_8patches.txt';

% Type of boundary conditions for each side of the domain
problem_data.nmnn_sides = [];
problem_data.drchlt_sides = [1 2 3 4 5 6];
problem_data.weak_drchlt_sides = [];

% Physical parameters
problem_data.c_diff = @(x, y) ones(size(x));
% problem_data.grad_c_diff = @(x, y) cat (1, ...
% reshape (zeros(size(x)), [1, size(x)]), ...
% reshape (zeros(size(x)), [1, size(x)]));

% Source and boundary terms
problem_data.f = @bilaplacian_rhs_Lshaped;
% problem_data.g = @(x, y, ind) bilaplacian_Lshaped_g_nmnn_8patches (x,y,ind);
problem_data.h = @(x, y, ind) solution_bilaplacian_Lshaped (x, y);

% Exact solution (optional)
problem_data.uex = @(x, y) solution_bilaplacian_Lshaped (x, y);
problem_data.graduex = @(x, y) solution_bilaplacian_Lshaped_grad (x, y);
problem_data.hessuex = @(x, y) solution_bilaplacian_Lshaped_hessian (x, y);

% DISCRETIZATION PARAMETERS
clear method_data
method_data.degree = [4 4]; % Degree of the splines
method_data.regularity = [2 2]; % Regularity of the splines
method_data.nsub = [16 16]; % Number of subdivisions of the coarsest mesh, with respect to the mesh in geometry
method_data.nquad = method_data.degree+1; % Points for the Gaussian quadrature rule

[geometry, msh, space, u] = mp_solve_bilaplace_C1 (problem_data, method_data);

[err_h2,~,~,err_h2s] = sp_h2_error (space, msh, u, problem_data.uex, problem_data.graduex, problem_data.hessuex)

subplot(1,2,1)
npts = [51 51];
sp_plot_solution (u, space, geometry, npts); shading interp; title('Computed solution')
vtk_pts = {linspace(0,1,npts(1)), linspace(0,1,npts(2))};
subplot(1,2,2)
for iptc = 1:msh.npatch
F = reshape (geometry(iptc).map(vtk_pts), [2, npts]);
X = reshape (F(1,:), npts); Y = reshape (F(2,:), npts);
surf(X, Y, problem_data.uex(X,Y));
hold on; shading interp
end
title('Exact solution')
56 changes: 56 additions & 0 deletions geopdes/inst/examples/base/ex_bilaplacian_hyperboloid.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
% PHYSICAL DATA OF THE PROBLEM
clear problem_data
problem_data.geo_name = 'geo_hyperboloid_ASG1.txt';

% Type of boundary conditions for each side of the domain
problem_data.nmnn_sides = [];
problem_data.drchlt_sides = [1 2 3 4 5 6];
problem_data.weak_drchlt_sides = [];

% Physical parameters
problem_data.c_diff = @(x, y, z) ones(size(x));

% Source and boundary terms
C = 1;
p = -1;
cons = 200;
problem_data.f = @(x, y, z) 16.*cons.*exp(1).^((-1).*cons.*(x.^2+y.^2+(x.^2+(-1).*y.^2).^2)).* ...
(1+4.*x.^2+4.*y.^2).^(-5).*(cons.^3.*(1+4.*x.^2+4.*y.^2).^3.*(4.* ...
x.^6+(-4).*x.^4.*((-1)+y.^2)+(y+2.*y.^3).^2+x.^2.*(1+8.*y.^2+(-4) ...
.*y.^4)).^2+8.*(64.*x.^8+x.^6.*(32+(-896).*y.^2)+2.*x.^4.*(7+(-48) ...
.*y.^2+832.*y.^4)+y.^2.*((-1)+14.*y.^2+32.*y.^4+64.*y.^6)+(-1).* ...
x.^2.*(1+60.*y.^2+96.*y.^4+896.*y.^6))+(-4).*cons.^2.*(1+4.*x.^2+ ...
4.*y.^2).^2.*(72.*x.^10+20.*x.^8.*(5+2.*y.^2)+x.^6.*(50+272.*y.^2+ ...
(-112).*y.^4)+(y+2.*y.^3).^2.*(1+7.*y.^2+18.*y.^4)+x.^4.*(11+142.* ...
y.^2+280.*y.^4+(-112).*y.^6)+x.^2.*(1+26.*y.^2+142.*y.^4+272.* ...
y.^6+40.*y.^8))+2.*cons.*(1+1568.*x.^10+19.*y.^2+174.*y.^4+796.* ...
y.^6+1752.*y.^8+1568.*y.^10+8.*x.^8.*(219+404.*y.^2)+4.*x.^6.*( ...
199+1064.*y.^2+2896.*y.^4)+2.*x.^4.*(87+786.*y.^2+3720.*y.^4+ ...
5792.*y.^6)+x.^2.*(19+244.*y.^2+1572.*y.^4+4256.*y.^6+3232.*y.^8)));
problem_data.h = @(x, y, z, ind) exp(-cons*(x.^2 + y.^2 + (x.^2 - y.^2).^2));

% Exact solution
problem_data.uex = @(x, y, z) exp(-cons*(x.^2 + y.^2 + (x.^2 - y.^2).^2));
problem_data.graduex = @(x, y, z) cat (1, ...
reshape ( (-2).*cons.*exp(1).^((-1).*cons.*(x.^2+y.^2+(x.^2+(-1).*y.^2).^2)).*x.*(1+4.*x.^2+4.*y.^2).^(-1).*(1+2.*x.^2+6.*y.^2), [1, size(x)]), ...
reshape ( (-2).*cons.*exp(1).^((-1).*cons.*(x.^2+y.^2+(x.^2+(-1).*y.^2).^2)).*y.*(1+6.*x.^2+2.*y.^2).*(1+4.*x.^2+4.*y.^2).^(-1), [1, size(x)]), ...
reshape ( (-4).*cons.*exp (1).^((-1).*cons.*(x.^2+y.^2+(x.^2+(-1).*y.^2).^2)).*(1+4.*x.^2+4.*y.^2).^(-1).*(x.^2+2.*x.^4+(-1).*y.^2+(-2).*y.^4), [1, size(x)]));

problem_data.lapuex = @(x, y, z) 4.*cons.*exp(1).^((-1).*cons.*(x.^2+y.^2+(x.^2+(-1).*y.^2).^2)).*( ...
1+4.*x.^2+4.*y.^2).^(-2).*((-1)+20.*cons.*x.^6+16.*cons.*x.^8+(( ...
-8)+cons).*y.^2+4.*((-5)+2.*cons).*y.^4+20.*cons.*y.^6+16.*cons.* ...
y.^8+x.^2.*((-8)+cons+(-24).*y.^2+16.*cons.*y.^2+44.*cons.*y.^4)+( ...
-4).*x.^4.*(5+cons.*((-2)+(-11).*y.^2+8.*y.^4)));

% DISCRETIZATION PARAMETERS
clear method_data
deg = 3;
method_data.degree = [deg deg]; % Degree of the splines
method_data.regularity = [deg-2 deg-2]; % Regularity of the splines
method_data.nsub = [16 16]; % Number of subdivisions of the coarsest mesh, with respect to the mesh in geometry
method_data.nquad = [deg+1 deg+1]; % Points for the Gaussian quadrature rule

[geometry, msh, space, u] = mp_solve_bilaplace_C1 (problem_data, method_data);
npts = [41 41];
sp_plot_solution (u, space, geometry); shading interp; axis equal
fprintf('Error in (equivalent) H2 seminorm = %g\n', sp_h2_equiv_lap_error(space, msh, u, problem_data.uex, problem_data.graduex, problem_data.lapuex));
Loading

0 comments on commit a233bc4

Please sign in to comment.