Commit 8b619052 authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Separate scalar- and vector-valued ContraintEdge value computation.

parent e5143d2f
......@@ -23,7 +23,7 @@ classdef ConstraintEdge < mafe.Constraint
if length(node) < 2
disp('ConstraintEdge: Error with given number of nodes < 2 !');
end
if length(cval) == 1 && length(type) == 1 % single scalar given
if length(cval) == 1 && length(obj.cdof) == 1 % single scalar given
obj.cval = cval * ones(1,length(node));
end
if length(cval) == 2 && length(obj.cdof) == 2 % single vector given
......@@ -76,19 +76,36 @@ classdef ConstraintEdge < mafe.Constraint
if ( L < eps )
continue;
end
% get parameters (in local coordinate system)
pT = self.cval( [ii, ii+1] );
pN = self.cval( [ii, ii+1]+length(self.node) );
% local element vector
vec_local = L/6.*[ +2, +1, 0, 0; ...
+1, +2, 0, 0; ...
0, 0, +2, +1; ...
0, 0, +1, +2 ] * [pT, pN]';
% global vector
vec_global = T' * vec_local;
%
idx = [ii, ii+1, ii+length(self.node), ii+1+length(self.node)];
vec( idx ) = vec( idx ) + vec_global';
switch (length(self.cdof))
% case: scalar-valued constraint
case 1
% get parameters (in local coordinate system)
p = self.cval( [ii, ii+1] );
% local element vector
vec_global = L/6.*[ +2, +1; ...
+1, +2] * [p]';
%
idx = [ii, ii+1];
vec( idx ) = vec( idx ) + vec_global';
% case: vector-valued constraint
case 2
% get parameters (in local coordinate system)
pT = self.cval( [ii, ii+1] );
pN = self.cval( [ii, ii+1]+length(self.node) );
% local element vector
vec_local = L/6.*[ +2, +1, 0, 0; ...
+1, +2, 0, 0; ...
0, 0, +2, +1; ...
0, 0, +1, +2 ] * [pT, pN]';
% global vector
vec_global = T' * vec_local;
%
idx = [ii, ii+1, ii+length(self.node), ii+1+length(self.node)];
vec( idx ) = vec( idx ) + vec_global';
% handle
otherwise
disp('ConstraintEdge: value() -- Error handling cases.')
end
end
end
end
......
......@@ -66,10 +66,10 @@ nn = nodes( [gl1.getGridNodes.id] ); % get nodes on GeoLine 4
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
% top edge: vertical distributed force q [N/m]
nn = nodes( [gl3.getGridNodes.id, gl5.getGridNodes.id] ); % get nodes on GeoLine 3+5
const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Disp3, mafe.DofType.Disp3], mafe.ConstraintType.Neumann, [ 0.0, -1.e4] );
const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Disp3], mafe.ConstraintType.Neumann, [ -1.e4] );
% right edge: vertical distributed force q [N/m]
nn = nodes( [gl10.getGridNodes.id] ); % get nodes on GeoLine 10
const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Disp3, mafe.DofType.Disp3], mafe.ConstraintType.Neumann, [ 0.0, -1.e4] );
const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Disp3], mafe.ConstraintType.Neumann, [ -1.e4] );
% --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
......
%% Example Plate 2e
% example on how to use the integrated mesh generator for
% a rectangular plate with circular hole
% loaded by lineload on the ring while completely fixed on the outer boundary
addpath([ '..' filesep '..' filesep '..' filesep ]);
%% Definition of sections
s1 = mafe.Section2D();
s1.E = 2.0e11; % [N/m^2]
s1.nue = 0.0; % [-]
s1.t = 1.0e-2; % [m]
% -----------------------------------------------------------MESH GENERATION---
%% Describe the goemetry of the problem for the mesh generator
% function that controls grid node distribution along geo lines
nf1 = mgen.MgNormFunction( mgen.MgNormFunctionType.QuadraticRev );
% create all geo points
gp1 = mgen.MgGeoPoint( [0.0, 0.0] );
gp2 = mgen.MgGeoPoint( [1.0, 0.0] );
gp3 = mgen.MgGeoPoint( [1.0, 1.0] );
gp4 = mgen.MgGeoPoint( [0.0, 1.0] );
% create all geo lines
gl1 = mgen.MgGeoLine2Point( [gp2, gp1], 21 ); % ask for 4 grid nodes along this line
gl2 = mgen.MgGeoLine2Point( [gp3, gp2], 21 ); % ask for 6 grid nodes along this line
gl3 = mgen.MgGeoLine2Point( [gp4, gp3], 2 ); % this will adapt automatically
gl4 = mgen.MgGeoLine2Point( [gp1, gp4], 2 ); % this will adapt automatically
% create all geo areas
ga1 = mgen.MgGeoArea4LineEllipsoidalVoid( [gl1, gl2, gl3, gl4], 20, 0.2, 0.2, 0.0, nf1 );
% create geo object and generate the mesh
geo = mgen.MgGeo( [gp1, gp2, gp3, gp4], [gl1, gl2, gl3, gl4], [ga1] );
% print geo and mesh data info
%fprintf(geo.print());
% plot mesh
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
geo.plot(); disp('Showing generated mesh. Press key to continue...'); pause
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of nodes in the system
nodes = mafe.Node.empty(0,0);
for gn = geo.getAllGridNodes()
nodes(end+1) = mafe.Node2D( gn.coord );
end
%% Definition of elements in the system
elems = mafe.Element.empty(0,0);
for ge = geo.getAllGridElems()
elems(end+1) = mafe.Plate2D( nodes( [ ge.nodes.id ] ), s1 );
end
%% Definition of constraints in the system
const = mafe.Constraint.empty(0,0);
% bottom edge: fix in x1 and x2
nn = nodes( [gl1.getGridNodes.id] ); % get nodes on GeoLine 1 (bottom)
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
% top edge: fix in x1 and x2
nn = nodes( [gl3.getGridNodes.id] ); % get nodes on GeoLine 3 (top)
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
% right edge: fix in x1 and x2
nn = nodes( [gl2.getGridNodes.id] ); % get nodes on GeoLine 2 (right)
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
% left edge: fix in x1 and x2
nn = nodes( [gl4.getGridNodes.id] ); % get nodes on GeoLine 4 (left)
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
% inner ring: apply line load
inner_ring_node_ids = [ga1.getGridNodesInner.id]; % get nodes on inner ring of GeoArea 1
nn = nodes( [inner_ring_node_ids, inner_ring_node_ids(1)] ); % make sure to append first node of closed line
const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Disp3], mafe.ConstraintType.Neumann, [-1.0e3] );
% --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
% %% Select an analysis type: static
ana = mafe.FeAnalysisStatic(fep);
% %% Calculate response
ana.analyse(); disp( sprintf('---> number dof = %5d', fep.ndofs) );
% -----------------------------------------------------------INSPECT RESULTS---
% %% Visualise system response
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
fep.plotSystem('reference');
fep.plotSystem('deformed', 1e2); view(45.0, 25.0);
fig = figure(2); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
fep.plotSystem('vonmises'); colorbar();
% fep.printSystem();
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment