Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Replies: 3   Last Post: Jan 31, 2011 11:00 AM

 Messages: [ Previous | Next ]
 Kwesi Moore Posts: 19 Registered: 11/27/09
Posted: Jan 7, 2011 1:16 PM

I have a quadrilateral mesh.i want to create refinements within the mesh.
This is what i have done so far....

element=[1 2 4 5;2 3 5 6;4 5 7 8;5 6 8 9];
coordinates=[-1 0 1 -1 0 1 -1 0 1;-1 -1 -1 0 0 0 1 1 1]';
dirichlet=[]1 2 1 4;2 3 4 7]';
neumann=[3 6 7 8;6 9 8 9]';
nodes2edge=[1 2 1 4 3 2 5 4 5 7 6 8;4 5 2 5 6 3 6 7 8 8 9 9]';

function [coordinates,NewElement,dirichlet,neumann,ParentElement]=...
refinement(coordinates, elements, dirichlet,neumann,nodes2edge)

% Compute number of edges, nodes and elements
no_of_edges = full(max(max(nodes2edge)));
no_of_nodes = size(coordinates,1);
no_of_elements = size(elements,1);

% Create New Arrays
Node2Edge = zeros(no_of_edges,2);
% Get the indices of the nodes
[a,b,c] = find(nodes2edge);
Node2Edge(c,:) = [a,b];

% Coordinates
coordinates(no_of_nodes+1:no_of_nodes+no_of_edges,:)=...
(coordinates(Node2Edge(:,1),:)+coordinates(Node2Edge(:,2),:))/2;

% Elements
NewElement = zeros(4*size(elements,1),3);
ParentElement = zeros(size(NewElement,1),1);

% Element Refinements
ell = 0;

for j = 1:no_of_elements
% describes the order of the faces >>left (face 1) right south north.
I = [nodes2edge(elements(j,1),elements(j,3));
nodes2edge(elements(j,2),elements(j,4));
nodes2edge(elements(j,1),elements(j,2));
nodes2edge(elements(j,3),elements(j,4))];
NewNode=(I+no_of_nodes)';
NewElement(ell+1:ell+5,:)=[NewNode([1 2 3 4]);
elements(j,1), NewNode(1), NewNode(2), NewNode(3);
NewNode(2), elements(j,2), NewNode(3), NewNode(4);
NewNode(3), NewNode(4), elements(j,3),NewNode(1);
NewNode(4), NewNode(1), NewNode(2), elements(j,4)];

ParentElement(ell+1:ell+5)=j;
ell=ell+5;
end

% Dirichlet and Neumann
tmp=0;
dirichletnr=diag(nodes2edge(dirichlet(:,1),dirichlet(:,2)));
for k=1:size(dirichlet,1)
dirichlet=[dirichlet(1:k-1+tmp,:);
dirichlet(k+tmp,1),dirichletnr(k)+no_of_nodes;
dirichletnr(k)+no_of_nodes,dirichlet(k+tmp,2);
dirichlet(k+1+tmp:size(dirichlet,1),:)];
tmp=tmp+1;
end

tmp=0;
if ~isempty(neumann)
neumannnr=diag(nodes2edge(neumann(:,1),neumann(:,2)));
end
for k=1:size(neumann,1)
neumann = [neumann(1:k-1+tmp,:);
neumann(k+tmp,1),neumannnr(k)+no_of_nodes;
neumannnr(k)+no_of_nodes,neumann(k+tmp,2);
neumann(k+1+tmp:size(neumann,1),:)];
tmp=tmp+1;
end

Date Subject Author
1/7/11 Kwesi Moore
1/7/11 Kwesi Moore
1/8/11 John D'Errico
1/31/11 Qwuasi Moore