Homec4science
Diffusion contact_internodes (master) Browse

rCONTACTINTERNODES/Code/Westergaard/setBoundaryConditions.m
fccc43844265master

Last Change
Raw File

setBoundaryConditions.m
View Options

function[Data]=setBoundaryConditions(Mesh, Data, varargin)
% setBoundaryConditions: Function which initializes the boundary data
% INPUT:
% Mesh: Structure containing all the mesh parameters
% Data: Structure containing all the data parameters
% OPTIONAL INPUT:
% Any type of boundary condition. If no Dirichlet boundary conditions are
% specified, the program returns an error
% Types of boundary conditions:
% Dirichlet
% Neumann
% OUTPUT:
% Data: Updated data structure with boundary conditions
%% Set boundary conditions
if mod(length(varargin{1}), 2)~=0
error('Missing a name or value')
end
% Retrieve boundary nodes
BC_tags=unique(Mesh.BC_tag);
% Dirichlet type
D_types={'dirichlet'};
% Neumann types
N_types={'neumann'};
% Available types of boundary conditions
available_types=union(D_types, N_types);
% Count to detect duplicate boundaries
count=[];
% Count number of Dirichlet and Neumann boundaries specified
nb_D_BC=0;
nb_N_BC=0;
% Retrieve labels and values entered by the user
labels_in=varargin{1}(1:2:end);
values_in=varargin{1}(2:2:end);
% Initilization
Data.BC=[];
% Checking for boundary conditions
c1=strcmp(labels_in, 'BC');
if any(c1)
BC=values_in{c1};
nb_BC=length(BC);
Data.BC=cell(nb_BC,1);
for i=1:nb_BC
bc=BC{i};
l=length(bc);
if mod(l,2)~=0
error('Missing a name or value')
end
for j=1:2:l
arg=lower(bc{j});
val=bc{j+1};
switch arg
case 'type'
if any(ismember(available_types, lower(val)))
Data.BC{i}.(arg)=lower(val);
if any(ismember(D_types, lower(val)))
nb_D_BC=nb_D_BC+1;
else
nb_N_BC=nb_N_BC+1;
end
else
error('Unrecognized boundary condition')
end
case 'edges'
edges=intersect(BC_tags, val);
if isempty(edges)
error('One of the edges specified does not exist')
elseif ~isempty(intersect(count, val)) || length(unique(val))~=length(val)
error('Duplicate edges detected')
end
count=[count val];
Data.BC{i}.(arg)=val;
case {'function'}
if isa(val,'double')
Data.BC{i}.label='const';
Data.BC{i}.(arg)= @(x) val;
elseif isa(val,'function_handle')
Data.BC{i}.label=detectDependency(val);
Data.BC{i}.(arg)=val;
else
error('Input parameters must be either constants or function handles')
end
case 'axis'
if strcmp(val, 'x') || strcmp(val, 'y') || strcmp(val, 'all')
Data.BC{i}.(arg)=val;
else
error('Unrecognized axis')
end
otherwise
error('Unrecognized argument')
end
end
end
end
%% Post-Processing: Partitioning the structure BC into Dirichlet BC and Neumann BC and checking for missing data
if nb_D_BC==0
error('Dirichlet boundary conditions must be prescribed')
end
BC=Data.BC;
l=length(BC);
% Initialize cells
D_BC=cell(nb_D_BC,1);
N_BC=cell(nb_N_BC,1);
% Set of Dirichlet and Neumann tags
D_tags=[];
N_tags=[];
% Initialize counters
d=1;
n=1;
for i=1:l
bc=BC{i};
type=bc.type;
if ~isfield(bc, 'edges')
error('The edges must be specified')
end
if any(ismember(D_types, type))
D_tags=[D_tags bc.edges];
D_BC{d}.type=type;
D_BC{d}.edges=bc.edges;
if isfield(bc, 'function')
D_BC{d}.function=bc.function;
D_BC{d}.label=bc.label;
else
warning('A function has not been prescribed and will be assumed zero')
D_BC{d}.function=@(x) [0; 0];
D_BC{d}.label='const';
end
if isfield(bc, 'axis')
D_BC{d}.axis=bc.axis;
else
D_BC{d}.axis='all';
end
d=d+1;
elseif any(ismember(N_types, type))
N_tags=[N_tags bc.edges];
N_BC{n}.type=type;
N_BC{n}.edges=bc.edges;
if isfield(bc, 'function')
N_BC{n}.function=bc.function;
N_BC{n}.label=bc.label;
else
warning('A function has not been prescribed and will be assumed zero')
N_BC{n}.function=@(x) [0; 0];
N_BC{n}.label='const';
end
if isfield(bc, 'axis')
N_BC{n}.axis=bc.axis;
else
N_BC{n}.axis='all';
end
n=n+1;
end
end
Data.D_BC=D_BC;
Data.N_BC=N_BC;
Data.D_tags=D_tags;
Data.N_tags=N_tags;
%% Sorting nodes
% If there is a node belonging to conflicting boundary conditions, it is
% set using the following priority list: Neumann, Dirichlet, Interior
BC_nodes=Mesh.BC_nodes;
BC_tag=Mesh.BC_tag;
nb_nodes=Mesh.nb_nodes;
Dof_set=Mesh.Dof_set;
set_D_nodes=cell(nb_D_BC,1);
set_N_nodes=cell(nb_N_BC,1);
set_Nd=cell(nb_D_BC,1);
set_Nn=cell(nb_N_BC,1);
for k=1:nb_D_BC
nodes=BC_nodes(BC_tag==D_BC{k}.edges,:);
set_D_nodes{k}=unique(nodes);
if strcmp(D_BC{k}.axis, 'x')
set_Nd{k}=Dof_set(1,set_D_nodes{k})';
elseif strcmp(D_BC{k}.axis, 'y')
set_Nd{k}=Dof_set(2,set_D_nodes{k})';
else
dof=Dof_set(:,set_D_nodes{k});
set_Nd{k}=dof(:);
end
end
for k=1:nb_N_BC
nodes=BC_nodes(BC_tag==N_BC{k}.edges,:);
set_N_nodes{k}=unique(nodes);
if strcmp(N_BC{k}.axis, 'x')
set_Nn{k}=Dof_set(1,set_N_nodes{k})';
elseif strcmp(N_BC{k}.axis, 'y')
set_Nn{k}=Dof_set(2,set_N_nodes{k})';
else
dof=Dof_set(:,set_N_nodes{k});
set_Nn{k}=dof(:);
end
end
% Define the sets of nodes
bc_nodes=unique(BC_nodes);
% Neumann nodes
Nodes_N=unique(cell2mat(set_N_nodes));
% Dirichlet nodes
Nodes_D=unique(cell2mat(set_D_nodes));
Nodes_D=setdiff(Nodes_D, Nodes_N);
Nodes_N=setdiff(bc_nodes, Nodes_D);
% Interior nodes
Nodes_I=setdiff(1:nb_nodes, union(Nodes_D, Nodes_N));
% Free nodes
Nodes_F=union(Nodes_I, Nodes_N);
% Define the sets of degrees of freedom
% Neumann nodes
Nn=unique(cell2mat(set_Nn));
% Dirichlet nodes
Nd=unique(cell2mat(set_Nd));
Nd=setdiff(Nd, Nn);
dof_n=Dof_set(:,bc_nodes);
dof_n=dof_n(:);
Nn=setdiff(dof_n, Nd);
% Interior nodes
Ni=setdiff(Dof_set(:), union(Nd, Nn));
% Free nodes
Nf=union(Ni, Nn);
% Correct the set for Dirichlet nodes
for k=1:nb_D_BC
set_D_nodes{k}=setdiff(set_D_nodes{k}, Nodes_N);
set_Nd{k}=setdiff(set_Nd{k}, Nn);
end
% Update data
Data.Nodes_I=Nodes_I;
Data.Nodes_N=Nodes_N;
Data.Nodes_D=Nodes_D;
Data.Nodes_F=Nodes_F;
Data.Ni=Ni;
Data.Nd=Nd;
Data.Nn=Nn;
Data.Nf=Nf;
Data.set_D_nodes=set_D_nodes;
Data.set_N_nodes=set_N_nodes;
Data.set_Nd=set_Nd;
Data.set_Nn=set_Nn;
%% Auxiliary function
function[label]=detectDependency(func)
string=func2str(func);
if ~strcmp(string(1:4), '@(x)')
error('Function handle must contain space variables')
end
string=string(5:end);
if contains(string, 'x') && contains(string, 't')
label='space_time';
elseif contains(string, 'x')
label='space';
elseif contains(string, 't')
label='time';
elseif contains(string(1), '0')
label='zero';
else
label='const';
end
end
end
c4science · Help