From a2769c04178de22f7bda7cb9884d4c1494149cbe Mon Sep 17 00:00:00 2001 From: "g.mandl" Date: Wed, 14 Jan 2026 18:42:17 +0100 Subject: [PATCH] Uploading ex6 --- Sheet_6/adapt_h.m | 18 ++++++++++++++++++ Sheet_6/adapt_r.m | 24 ++++++++++++++++++++++++ Sheet_6/assembling.m | 20 ++++++++++++++++++++ Sheet_6/jumps_flux.m | 16 ++++++++++++++++ 4 files changed, 78 insertions(+) create mode 100644 Sheet_6/adapt_h.m create mode 100644 Sheet_6/adapt_r.m create mode 100644 Sheet_6/assembling.m create mode 100644 Sheet_6/jumps_flux.m diff --git a/Sheet_6/adapt_h.m b/Sheet_6/adapt_h.m new file mode 100644 index 0000000..1e4e363 --- /dev/null +++ b/Sheet_6/adapt_h.m @@ -0,0 +1,18 @@ +% h-adaptivity +% for a given mesh and solution we generate a new adapted mesh by splitting +% the elements with large errors into two new elements + +function mesh = adapt_h(nodes,u,lambda) + n_nodes = length(nodes); + flux_jumps = jumps_flux(nodes,u,lambda); + alpha = 0.5; % parameter for choosing elements to refine + crit_error = alpha*max(abs(flux_jumps)); + mesh = nodes; % will be the new nodes/mesh + % finding elements to refine + for i=2:n_nodes-1 + if abs(flux_jumps(i)) > crit_error + mesh = [mesh, nodes(i-1)+(nodes(i)-nodes(i-1))/2, nodes(i)+(nodes(i+1)-nodes(i))/2]; + end + end + mesh = unique(mesh); +end \ No newline at end of file diff --git a/Sheet_6/adapt_r.m b/Sheet_6/adapt_r.m new file mode 100644 index 0000000..85f49df --- /dev/null +++ b/Sheet_6/adapt_r.m @@ -0,0 +1,24 @@ +% r-adaptivity +% for a given mesh and solution we generate a new adapted mesh by moving +% the existing nodes within the mesh; they get moved to a position in order +% to equally distribute the error over the intervall +% we use the De Boor's algorithm (Huang, Russell; Adaptive Moving Mesh +% Methods; $ 2.2.1) + +function mesh = adapt_r(nodes,u,lambda) + n_nodes = length(nodes); + flux_jumps = abs(jumps_flux(nodes,u,lambda)); + p = 1/2*(flux_jumps(1:end-1) + flux_jumps(2:end)); % has values for each element + h_vec = nodes(2:end) - nodes(1:end-1); + P = cumsum(h_vec'.*p); + P = [0;P]; + xi = linspace(0,1,n_nodes); + mesh = nodes; + for j=2:n_nodes-1 + idx_k = find(xi(j)*P(end) <= P, 1, 'first'); + if idx_k > 1 + x_j = nodes(idx_k-1) + xi(j)*(P(end)-P(idx_k -1))/p(idx_k -1); + mesh(j) = x_j; + end + end +end \ No newline at end of file diff --git a/Sheet_6/assembling.m b/Sheet_6/assembling.m new file mode 100644 index 0000000..947d15f --- /dev/null +++ b/Sheet_6/assembling.m @@ -0,0 +1,20 @@ +% assembling of the stiffness matrix and the load vector for a given mesh +% in 1D + +function [K,f_vec] = assembling(nodes,lambda,f) + n_nodes = length(nodes); + K = zeros(n_nodes); + f_vec = zeros(n_nodes,1); + h_diff_vec = nodes(2:end) - nodes(1:end-1); % step width + for k = 1:n_nodes-1 + % stiffness matrix + lambda_int = integral(lambda, nodes(k), nodes(k+1)); + K_loc = lambda_int/h_diff_vec(k)^2*[1,-1;-1,1]; + K(k:k+1,k:k+1) = K(k:k+1,k:k+1) + K_loc; + % right hand side + phi_left = @(x) (nodes(k+1)-x)/h_diff_vec(k).*f(x); + phi_right = @(x) (x-nodes(k))/h_diff_vec(k).*f(x); + f_loc = [integral(phi_left,nodes(k),nodes(k+1)); integral(phi_right,nodes(k),nodes(k+1))]; + f_vec(k:k+1) = f_vec(k:k+1) + f_loc; + end +end \ No newline at end of file diff --git a/Sheet_6/jumps_flux.m b/Sheet_6/jumps_flux.m new file mode 100644 index 0000000..0a4b95b --- /dev/null +++ b/Sheet_6/jumps_flux.m @@ -0,0 +1,16 @@ +% calculation of the flux jumps according to the chosen mesh +% nodes - corresponding to the mesh +% u - approximation +% flux_jumps - jump of the flux in each node of the mesh including the +% paramter function lambda + +function flux_jumps = jumps_flux(nodes, u, lambda) + n_nodes = length(nodes); + flux_jumps = zeros(n_nodes,1); + diff_u = u(2:end) - u(1:end-1); + diff_x = nodes(2:end) - nodes(1:end-1); + eps = 1e-8; + for i = 2:n_nodes-1 + flux_jumps(i) = diff_u(i)/diff_x(i)*lambda(nodes(i)+eps) - diff_u(i-1)/diff_x(i-1)*lambda(nodes(i)-eps); + end +end \ No newline at end of file