Skip to content

Commit

Permalink
New functions, intead of int_boundary_term
Browse files Browse the repository at this point in the history
  • Loading branch information
@rafavzqz
rafavzqz committed Jul 22, 2024
1 parent c4960b7 commit 32e9309
Show file tree
Hide file tree
Showing 5 changed files with 226 additions and 98 deletions.
46 changes: 25 additions & 21 deletions geopdes/INDEX
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,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 @@ -183,27 +185,29 @@ Multipatch: functions of different kinds to solve problems in domains defined by
@sp_multipatch/sp_weak_drchlt_bc_stokes
@sp_multipatch/sp_refine
sp_multipatch_C1
@sp_multipatch/sp_l2_error
@sp_multipatch/sp_h1_error
@sp_multipatch/sp_h2_error
@sp_multipatch/sp_h2_equiv_lap_error
@sp_multipatch/sp_plot_solution
@sp_multipatch/sp_to_vtk
@sp_multipatch/sp_eval_phys
@sp_multipatch/sp_evaluate_element_list
@sp_multipatch/sp_get_cells
@sp_multipatch/sp_get_basis_functions
@sp_multipatch/sp_get_functions_on_patch
@sp_multipatch/sp_get_local_interior_functions
@sp_multipatch/sp_get_neighbors
@sp_multipatch/sp_compute_Cpatch
@sp_multipatch/sp_compute_Cpatch_vector
@sp_multipatch/sp_bilaplacian_drchlt_C1_exact
@sp_multipatch/sp_bilaplacian_drchlt_C1
@sp_multipatch/sp_drchlt_C1_shells
@sp_multipatch/sp_nitsche_KL_rotation
@sp_multipatch/sp_weak_drchlt_bc_laplace
@sp_multipatch/sp_refine
@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
Expand Down
68 changes: 68 additions & 0 deletions geopdes/inst/multipatch/@sp_multipatch_C1/op_nitsche_consistency_cahn_hilliard.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
% OP_NITSCHE_CONSISTENCY_CAHN_HILLIARD: assemble the matrix A = [a(i,j)], a(i,j) = (epsilon (grad v n)_j, Delta u_i),
% with n the normal vector to the boundary, using the same space for trial and test functions.
%
% mat = op_nitsche_consistency_cahn_hilliard (space, msh, sides, [coeff]);
%
% INPUT:
%
% space: object representing the space of trial functions (see sp_multipatch_C1)
% msh: object defining the mesh (see msh_multipatch)
% sides: boundary sides on which to compute the integrals
% coeff: function handle to compute the epsilon coefficient. If empty, it is taken equal to one.
%
% OUTPUT:
%
% mat: assembled matrix
%
% Copyright (C) 2023, 2024 Michele Torre, Rafael Vazquez
%
% 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
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.

% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.

function [A] = op_nitsche_consistency_cahn_hilliard (space, msh, nmnn_sides, coeff)

if (nargin == 3 || isempty(coeff))
coeff = @(varargin) ones (size(varargin{1}));
end

if (~isempty(nmnn_sides))

A = spalloc (space.ndof, space.ndof, 3*space.ndof);

for iref = nmnn_sides
for bnd_side = 1:msh.boundaries(iref).nsides
iptc = msh.boundaries(iref).patches(bnd_side);
iside = msh.boundaries(iref).faces(bnd_side);

msh_side = msh_eval_boundary_side (msh.msh_patch{iptc}, iside ) ;
msh_side_int = msh_boundary_side_from_interior (msh.msh_patch{iptc}, iside ) ;

sp_side = space.sp_patch{iptc}.constructor (msh_side_int);
sp_side = sp_precompute ( sp_side , msh_side_int , 'gradient' , true, 'laplacian', true );

for idim = 1:msh.rdim
x{idim} = reshape (msh_side.geo_map(idim,:,:), msh_side.nqn, msh_side.nel);
end
coe_side = coeff (x{:});
tmp = op_gradv_n_laplaceu(sp_side ,sp_side ,msh_side, coe_side);

[Cpatch, Cpatch_cols] = sp_compute_Cpatch (space, iptc);
A(Cpatch_cols,Cpatch_cols) = ...
A(Cpatch_cols,Cpatch_cols) + Cpatch.' * tmp * Cpatch;
end
end

else
A = [];
end
end
84 changes: 41 additions & 43 deletions geopdes/inst/multipatch/mp_solve_cahn_hilliard_C1.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,6 @@
%
% problem_data: a structure with data of the problem. It contains the fields:
% - geo_name: name of the file containing the geometry
% - periodic_directions: parametric directions along which to apply periodic conditions (may be empty)
% - lambda: parameter representing the length scale of the problem, and the width of the interface
% - mu: function handle to compute mu (from the double well function)
% - dmu: function handle to compute the derivative of mu
Expand Down Expand Up @@ -50,7 +49,7 @@
% Only periodic and Neumann boundary conditions are implemented. Neumann
% conditions are considered by default.
%
% Copyright (C) 2023 Michele Torre, Rafael Vazquez
% Copyright (C) 2023, 2024 Michele Torre, Rafael Vazquez
%
% 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 Down Expand Up @@ -125,12 +124,11 @@
lapl_mat = op_laplaceu_laplacev_mp (space, space, msh, lambda);

% Compute the boundary term
bnd_mat = int_boundary_term (space, msh, lambda, nmnn_sides);
bnd_mat = op_nitsche_consistency_cahn_hilliard (space, msh, nmnn_sides, lambda);

% Compute the penalty matrix
[Pen, pen_rhs] = op_penalty_dudn (space, msh, nmnn_sides, Cpen_nitsche);


%%-------------------------------------------------------------------------
% Initial conditions
time = 0;
Expand Down Expand Up @@ -292,51 +290,51 @@
% Tangent stiffness matrix (mass is not considered here)
stiff_mat = term2 + term2K + lapl_mat;

% in case of neumann BC, add boundary terms
% In case of Neumann BC, add boundary terms
if (~isempty(bnd_mat))
Res_gl = Res_gl - (bnd_mat + bnd_mat.') * u_a + Pen*u_a - pen_rhs;
stiff_mat = stiff_mat - (bnd_mat + bnd_mat.') + Pen;
end
end

%--------------------------------------------------------------------------
% Boundary term, \int_\Gamma (\Delta u) (\partial v / \partial n)
%--------------------------------------------------------------------------

function [A] = int_boundary_term (space, msh, lambda, nmnn_sides)

if (~isempty(nmnn_sides))

A = spalloc (space.ndof, space.ndof, 3*space.ndof);

for iref = nmnn_sides
for bnd_side = 1:msh.boundaries(iref).nsides
iptc = msh.boundaries(iref).patches(bnd_side);
iside = msh.boundaries(iref).faces(bnd_side);

msh_side = msh_eval_boundary_side (msh.msh_patch{iptc}, iside ) ;
msh_side_int = msh_boundary_side_from_interior (msh.msh_patch{iptc}, iside ) ;

sp_side = space.sp_patch{iptc}.constructor (msh_side_int);
sp_side = sp_precompute ( sp_side , msh_side_int , 'gradient' , true, 'laplacian', true );

for idim = 1:msh.rdim
x{idim} = reshape (msh_side.geo_map(idim,:,:), msh_side.nqn, msh_side.nel);
end
coe_side = lambda (x{:});
tmp = op_gradv_n_laplaceu(sp_side ,sp_side ,msh_side, coe_side);

[Cpatch, Cpatch_cols] = sp_compute_Cpatch (space, iptc);
A(Cpatch_cols,Cpatch_cols) = ...
A(Cpatch_cols,Cpatch_cols) + Cpatch.' * tmp * Cpatch;
end
end

else
A = [];
end
end

% %--------------------------------------------------------------------------
% % Boundary term, \int_\Gamma (\Delta u) (\partial v / \partial n)
% %--------------------------------------------------------------------------
%
% function [A] = int_boundary_term (space, msh, lambda, nmnn_sides)
%
% if (~isempty(nmnn_sides))
%
% A = spalloc (space.ndof, space.ndof, 3*space.ndof);
%
% for iref = nmnn_sides
% for bnd_side = 1:msh.boundaries(iref).nsides
% iptc = msh.boundaries(iref).patches(bnd_side);
% iside = msh.boundaries(iref).faces(bnd_side);
%
% msh_side = msh_eval_boundary_side (msh.msh_patch{iptc}, iside ) ;
% msh_side_int = msh_boundary_side_from_interior (msh.msh_patch{iptc}, iside ) ;
%
% sp_side = space.sp_patch{iptc}.constructor (msh_side_int);
% sp_side = sp_precompute ( sp_side , msh_side_int , 'gradient' , true, 'laplacian', true );
%
% for idim = 1:msh.rdim
% x{idim} = reshape (msh_side.geo_map(idim,:,:), msh_side.nqn, msh_side.nel);
% end
% coe_side = lambda (x{:});
% tmp = op_gradv_n_laplaceu(sp_side ,sp_side ,msh_side, coe_side);
%
% [Cpatch, Cpatch_cols] = sp_compute_Cpatch (space, iptc);
% A(Cpatch_cols,Cpatch_cols) = ...
% A(Cpatch_cols,Cpatch_cols) + Cpatch.' * tmp * Cpatch;
% end
% end
%
% else
% A = [];
% end
% end
%
% %--------------------------------------------------------------------------
% % penalty term
% %--------------------------------------------------------------------------
Expand Down
66 changes: 32 additions & 34 deletions geopdes/inst/solve/solve_cahn_hilliard.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@
% Only periodic and Neumann boundary conditions are implemented. Neumann
% conditions are considered by default.
%
% Copyright (C) 2023 Michele Torre, Rafael Vazquez
% Copyright (C) 2023, 2024 Michele Torre, Rafael Vazquez
%
% 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 Down Expand Up @@ -133,13 +133,12 @@
lapl_mat = op_laplaceu_laplacev_tp (space, space, msh, lambda);

% Compute the boundary term
bnd_mat = int_boundary_term (space, msh, lambda, nmnn_sides);
bnd_mat = op_nitsche_consistency_cahn_hilliard (space, msh, nmnn_sides, lambda);

% Compute the penalty matrix
%[Pen, pen_rhs] = penalty_matrix (space, msh, nmnn_sides, Cpen_nitsche);
[Pen, pen_rhs] = op_penalty_dudn (space, msh, nmnn_sides, Cpen_nitsche);


%%-------------------------------------------------------------------------
% Initial conditions
time = 0;
Expand Down Expand Up @@ -260,7 +259,6 @@
end
norm_res = norm(Res_gl);


if (norm_res/norm_res_0 < tol_rel_res)
disp(strcat('iteration n°=',num2str(iter)))
disp(strcat('norm (abs) residual=',num2str(norm_res)))
Expand Down Expand Up @@ -302,42 +300,42 @@
% Tangent stiffness matrix (mass is not considered here)
stiff_mat = term2 + term2K + lapl_mat;

% In case of neumann BC, add boundary terms
% In case of Neumann BC, add boundary terms
if (~isempty(bnd_mat))
Res_gl = Res_gl - (bnd_mat + bnd_mat.') * u_a + Pen*u_a - pen_rhs;
stiff_mat = stiff_mat - (bnd_mat + bnd_mat.') + Pen;
end
end

%--------------------------------------------------------------------------
% Boundary term, \int_\Gamma (\Delta u) (\partial v / \partial n)
%--------------------------------------------------------------------------

function [A] = int_boundary_term (space, msh, lambda, nmnn_sides)

if (~isempty(nmnn_sides))

A = spalloc (space.ndof, space.ndof, 3*space.ndof);

for iside = 1:numel(nmnn_sides)

msh_side = msh_eval_boundary_side (msh, nmnn_sides(iside) ) ;
msh_side_int = msh_boundary_side_from_interior (msh, nmnn_sides(iside) ) ;
sp_side = space.constructor ( msh_side_int) ;
sp_side = sp_precompute ( sp_side , msh_side_int , 'gradient' , true, 'laplacian', true );

for idim = 1:msh.rdim
x{idim} = reshape (msh_side.geo_map(idim,:,:), msh_side.nqn, msh_side.nel);
end
coe_side = lambda (x{:});

A = A + op_gradv_n_laplaceu(sp_side ,sp_side ,msh_side, coe_side);
end

else
A = [];
end
end
% %--------------------------------------------------------------------------
% % Boundary term, \int_\Gamma (\Delta u) (\partial v / \partial n)
% %--------------------------------------------------------------------------
%
% function [A] = int_boundary_term (space, msh, lambda, nmnn_sides)
%
% if (~isempty(nmnn_sides))
%
% A = spalloc (space.ndof, space.ndof, 3*space.ndof);
%
% for iside = 1:numel(nmnn_sides)
%
% msh_side = msh_eval_boundary_side (msh, nmnn_sides(iside) ) ;
% msh_side_int = msh_boundary_side_from_interior (msh, nmnn_sides(iside) ) ;
% sp_side = space.constructor ( msh_side_int) ;
% sp_side = sp_precompute ( sp_side , msh_side_int , 'gradient' , true, 'laplacian', true );
%
% for idim = 1:msh.rdim
% x{idim} = reshape (msh_side.geo_map(idim,:,:), msh_side.nqn, msh_side.nel);
% end
% coe_side = lambda (x{:});
%
% A = A + op_gradv_n_laplaceu(sp_side ,sp_side ,msh_side, coe_side);
% end
%
% else
% A = [];
% end
% end

% %--------------------------------------------------------------------------
% % Penalty term
Expand Down
Loading

0 comments on commit 32e9309

Please sign in to comment.