Commit e5143d2f authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Add access to grid nodes and elements to MgGeoArea.

parent 9cab594e
classdef MgGeoArea < mgen.MgGeoObject
% An (unspecified) area of the geometry
properties
properties (Access = protected)
gridNodes = mgen.MgGridNode.empty(0,0);
gridElems = mgen.MgGridElem.empty(0,0);
end
methods (Abstract)
prepareLines(obj)
makeGridNodes(obj, gridnodes)
makeGridElems(obj, gridelems)
getGridNodes(obj)
end
methods
function gn = getGridNodes(self)
gn = self.gridNodes;
end
function ge = getGridElems(self)
ge = self.gridElems;
end
end
end
......@@ -3,9 +3,6 @@ classdef MgGeoArea4Line < mgen.MgGeoArea
properties (Access = protected)
geoLines = mgen.MgGeoLine.empty(0,4);
gridNodes = mgen.MgGridNode.empty(0,0);
gridElems = mgen.MgGridElem.empty(0,0);
end
methods
......@@ -67,7 +64,7 @@ classdef MgGeoArea4Line < mgen.MgGeoArea
self.geoLines(4).switchStartGridNode();
eps = 1.e-6;
% create grid nodes
for ii = 2:dim0-1 % straight edges
......@@ -75,17 +72,17 @@ classdef MgGeoArea4Line < mgen.MgGeoArea
A = self.geoLines(1).getGridNode( 1).coord;
B = self.geoLines(1).getGridNode(dim0).coord;
C = self.geoLines(1).getGridNode( ii).coord;
z1 = norm((C-A)/(B-A));
A = self.geoLines(3).getGridNode( 1).coord;
B = self.geoLines(3).getGridNode(dim0).coord;
C = self.geoLines(3).getGridNode( ii).coord;
z3 = norm((C-A)/(B-A));
pos = (ii-1)*dim1 + 2;
for jj = 2:dim1-1 % curved edges
% determine local coordinates of "same" points on both straight edges
......@@ -100,17 +97,17 @@ classdef MgGeoArea4Line < mgen.MgGeoArea
C = self.geoLines(4).getGridNode( jj).coord;
z4 = norm((C-A)/(B-A));
A = self.geoLines(2).getGridNode(jj).coord;
B = self.geoLines(4).getGridNode(jj).coord;
C = self.geoLines(1).getGridNode(ii).coord;
D = self.geoLines(3).getGridNode(ii).coord;
zi = (ii-1)/(dim0-1); % local z along curve
zj = (jj-1)/(dim1-1); % local z along curve
zzi = (1-zi) * z2 + zi * z4; % interpolate outer pos
zzj = (1-zj) * z1 + zj * z3;
zzj = (1-zj) * z1 + zj * z3;
N = ( (1-zzj)*A + (zzj)*B + (1-zzi)*C + (zzi)*D )/2;
% position node
......@@ -122,7 +119,7 @@ classdef MgGeoArea4Line < mgen.MgGeoArea
pos = pos + 1;
end
end
% sort in existing nodes
for ii = 1:dim1
self.gridNodes(ii) = self.geoLines(2).getGridNode(ii);
......@@ -169,9 +166,5 @@ classdef MgGeoArea4Line < mgen.MgGeoArea
end
end
function gn = getGridNodes(self)
gn = self.gridNodes;
end
end
end
classdef MgGeoAreaCircle < mgen.MgGeoArea
% A circular area of the (disk) geometry
properties
properties (Access = protected)
nrings = 2; % the number of rings between center part and outer line
ncircu = 2; % the number of rings between center part and outer line
recfac = 0.3; % fraction of the radius to be accupied by internal quadrilateral
center; % the center point of the circle
gridNodes = mgen.MgGridNode.empty(0,0);
gridElems = mgen.MgGridElem.empty(0,0);
ringNodes; % generated grid nodes on the outer boundary (closed connectivity)
end
......@@ -62,10 +59,6 @@ classdef MgGeoAreaCircle < mgen.MgGeoArea
% provide your code that generates the grid elements here and add them to gev (the container holding all grid elements)
end
function gn = getGridNodes(self)
gn = self.gridNodes;
end
function gn = getGridNodesOuter(self)
gn = self.ringNodes;
end
......
%% Example Membrane in Plane Stress 3a
% example on how to use the integrated mesh generator for
% a sblock under constant distributed vertical load q along the top edge
% composed of two materials
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]
s1.state = 'plane_stress';
s2 = mafe.Section2D();
s2.E = 1.0e11; % [N/m^2]
s2.nue = 0.3; % [-]
s2.t = 1.0e-2; % [m]
s2.state = 'plane_stress';
% -----------------------------------------------------------MESH GENERATION---
%% Describe the goemetry of the problem for the mesh generator
% create all geo points
gp1 = mgen.MgGeoPoint( [0.0, 0.0] );
gp2 = mgen.MgGeoPoint( [1.0, 0.0] );
gp3 = mgen.MgGeoPoint( [1.0, 0.5] );
gp4 = mgen.MgGeoPoint( [0.0, 0.5] );
gp5 = mgen.MgGeoPoint( [1.0, 1.0] );
gp6 = mgen.MgGeoPoint( [0.0, 1.0] );
% create all geo lines
gl1 = mgen.MgGeoLine2Point( [gp1, gp2], 8 ); % ask for 4 grid nodes along this line
gl2 = mgen.MgGeoLine2Point( [gp2, gp3], 4 ); % ask for 6 grid nodes along this line
gl3 = mgen.MgGeoLine2Point( [gp3, gp4], 2 ); % this will adapt automatically
gl4 = mgen.MgGeoLine2Point( [gp4, gp1], 2 ); % this will adapt automatically
gl5 = mgen.MgGeoLine2Point( [gp3, gp5], 16 ); % this will adapt automatically
gl6 = mgen.MgGeoLine2Point( [gp5, gp6], 2 ); % this will adapt automatically
gl7 = mgen.MgGeoLine2Point( [gp6, gp4], 2 ); % this will adapt automatically
% create all geo areas
ga1 = mgen.MgGeoArea4Line( [gl1, gl2, gl3, gl4] );
ga2 = mgen.MgGeoArea4Line( [gl3, gl5, gl6, gl7] );
% create geo object and generate the mesh
geo = mgen.MgGeo( [gp1, gp2, gp3, gp4, gp5, gp6], [gl1, gl2, gl3, gl4, gl5, gl6, gl7], [ga1, ga2] );
% 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 blocks
elems = mafe.Element.empty(0,0);
for ge = ga1.getGridElems() % define s1 for lower block
elems(end+1) = mafe.Membrane2D( nodes( [ ge.nodes.id ] ), s1 );
end
for ge = ga2.getGridElems() % define s2 for upper block
elems(end+1) = mafe.Membrane2D( nodes( [ ge.nodes.id ] ), s2 );
end
%% Definition of constraints in the system
const = mafe.Constraint.empty(0,0);
% bottom edge: fix in x2
nn = nodes( [gl1.getGridNodes.id] ); % get nodes on GeoLine 1 (bottom)
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% left bottom corner node: fix in x1
nn = nodes( [gp1.getGridNodes.id] ); % get node at GeoPoint 1 (bottom left)
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
% top edge: upwards distributed force q = 1e3 [N/m] translates to edge normal load
nn = nodes( [gl6.getGridNodes.id] ); % get nodes on GeoLine 3 (top)
const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Disp1, mafe.DofType.Disp2], mafe.ConstraintType.Neumann, [0.0, 1.e3] );
% --------------------------------------------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', 1e5);
%fep.plotSystem('stress22'); 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