stokes_compute_drag_and_lift.m

Go to the documentation of this file.

function [Cd, Cl] = stokes_compute_drag_and_lift(model, model_data, sim_data)
%function [Cd, Cl] = stokes_compute_drag_and_lift(model, model_data, sim_data)
% Drag and lift coefficients are computed via volume integrals.
%
% Required fields of model:
%  inner_boundary_rectangle: vector [xleft xright yleft yright]
%  representing a rectangle, which contains the boundary
%  U_mean: mean inflow velocity
%  L: characteristic length

% get dofs at inner boundary
dirichlet_dofs = model_data.bc_info.dirichlet_gids;

v = Fem.DiscFunc([], model_data.df_info);
global_coordinates = v.dofs_global_coordinate(dirichlet_dofs, :);

obstacle_boundary = global_coordinates(:, 1) > model.inner_boundary_rectangle(1) ...
    & global_coordinates(:, 1) < model.inner_boundary_rectangle(2) ...
    & global_coordinates(:, 2) > model.inner_boundary_rectangle(3) ...
    & global_coordinates(:, 2) < model.inner_boundary_rectangle(4);

dirichlet_dofs = dirichlet_dofs(obstacle_boundary);

% delete dofs for operator computation
old_dirichlet_gids = model_data.bc_info.dirichlet_gids;
model_data.bc_info.dirichlet_gids = setdiff(model_data.bc_info.dirichlet_gids, ...
    dirichlet_dofs);
model.has_dirichlet_values = ~isempty(model_data.bc_info.dirichlet_gids);
if model.has_dirichlet_values
    model.decomp_mode = 1;
    model_data.bc_info.dirichlet_dof_vector_components = ...
        Fem.Assembly.dirichlet_dof_vector(model, model_data.df_info, model_data.bc_info);
    model.decomp_mode = 0;
end

[A, f] = model.operators(model, model_data);
if model.has_nonlinearity
    model.uh = sim_data.uh;
    A = A + model.operators(model, model_data);
end
% undo changes
model_data.bc_info.dirichlet_gids = old_dirichlet_gids;
if model.has_dirichlet_values
    model.decomp_mode = 1;
    model_data.bc_info.dirichlet_dof_vector_components = ...
        Fem.Assembly.dirichlet_dof_vector(model, model_data.df_info, model_data.bc_info);
    model.decomp_mode = 0;
end

% drag
drag_dofs = intersect(dirichlet_dofs, get_dof_ids_global(model_data.df_info, [1 1]));
v.dofs(drag_dofs) = 1;

Cd = -2/model.U_mean(model)/model.U_mean(model)/model.L(model) * v.dofs'*(A*sim_data.uh.dofs - f);

% lift
lift_dofs = intersect(dirichlet_dofs, get_dof_ids_global(model_data.df_info, [1 2]));
v.dofs = zeros(v.ndofs, 1);
v.dofs(lift_dofs) = 1;

Cl = -2/model.U_mean(model)/model.U_mean(model)/model.L(model) * v.dofs'*(A*sim_data.uh.dofs - f);

end