3 % structure representing the fem-space information shared by all
4 % fem-functions. Implemented as handle
class, in order to be linked
7 % B. Haasdonk 13.1.2011
9 % adapted to
new assembly
17 % dimension of ranges space
20 % index vector of global DOFs
26 % number of DOFs per grid element
29 %
object of type .gridbase
32 % number of lagrange node indices
35 % the lagrange points on nodes
38 % the lagrange points on edges
41 % determinant of reference mapping
44 % global IDs of dirichlet DOFs
47 % element IDs of Robin DOFs
50 % edge IDs of Robin DOFs
53 % element IDs of Neumann DOFs
56 % edge IDs of Neumann DOFs
59 % weight matrix
for `L^2` inner product
60 l2_inner_product_matrix;
62 % weight matrix
for `H^1_0` inner product
63 h10_inner_product_matrix;
65 % weight matrix
for regularized `H^1_0` inner product
66 regularized_h10_inner_product_matrix;
73 % required fields of model:
74 % pdeg: polynomial degree
75 df_info.pdeg = model.pdeg;
76 df_info.dimrange = model.dimrange;
77 df_info.ndofs = fem_ndofs(model,grid);
78 df_info.ndofs_per_element = fem_ndofs_per_element(model);
82 df_info.lagrange_nodes = lagrange_nodes_lcoord(model.pdeg);
83 df_info.nlagrange_nodes = size(df_info.lagrange_nodes,1);
84 df_info.lagrange_edges = lagrange_nodes_edges(model.pdeg);
85 df_info.detDF = grid.A(:)*2; % determinant of reference mapping per element;
86 df_info.dirichlet_gids = [];
88 % l2_inner_product and h10_inner product
for now only
for
89 % scalar functions. Extend later.
90 % make sure, that the folliwing is called with empty
91 % dirichlet nodes, as fem_matrix_volume_component_assembly does
92 % dirichlet treatment!!!
93 % reaction term == 1 then reaction matrix is l2matrix
94 if df_info.dimrange == 1
95 model.reaction = @(grid,eindices,loc,params) ones(size(eindices,1),1);
96 model.qdeg = 2*df_info.pdeg;
97 df_info.l2_inner_product_matrix = ...
98 fem_matrix_volume_part_assembly(@fem_matrix_reac_integral_kernel,model,df_info);
99 % H10_inner_product_matrix: identity diffusion
100 model.diffusivity_tensor = @(grid,eindices,loc,params) ...
101 [ones(size(eindices,1),1), zeros(size(eindices,1),1), ...
102 zeros(size(eindices,1),1), ones(size(eindices,1),1)];
103 df_info.h10_inner_product_matrix = ...
104 fem_matrix_volume_part_assembly(@fem_matrix_diff_integral_kernel,model,df_info);
107 % find dirichlet indices as preparation
108 is_dirichlet_gid = zeros(df_info.ndofs,1);
109 % dirichlet-treatment: insert unit-vector row into A in dirichlet
110 % rows. Little loop over local lagrange node index
111 %
this is repeated in r assembly as these will be kept separated
112 for i=1:df_info.nlagrange_nodes
113 % search all elements, in which the i-th Lagrange node is
115 eids = find(df_info.lagrange_edges(i,:)==1);
117 for e = eids; % possibly one point belongs to 2 edges!!
118 elids = find(grid.NBI(:,e)==-1);
119 gids = df_info.global_dof_index(elids,i);
120 is_dirichlet_gid(gids) = 1;
124 df_info.dirichlet_gids = find(is_dirichlet_gid);
126 % find robin boundary elements and edge numbers
127 ind = find(grid.NBI==-3);
128 [df_info.robin_elind,df_info.robin_edgeind] = ind2sub(size(grid.NBI),ind);
130 % find neumann boundary elements and edge numbers
131 ind = find(grid.NBI==-2);
132 [df_info.neumann_elind,df_info.neumann_edgeind] = ind2sub(size(grid.NBI),ind);
133 % regularized_h10_inner_product_matrix: make unit rows in
136 if df_info.dimrange == 1
137 reg_mat = df_info.h10_inner_product_matrix;
138 if ~isempty(df_info.dirichlet_gids)
139 reg_mat(df_info.dirichlet_gids,:) = 0;
140 reg_mat(:,df_info.dirichlet_gids) = 0;
141 A_dirichlet = sparse(df_info.dirichlet_gids,...
142 df_info.dirichlet_gids, ...
143 ones(size(df_info.dirichlet_gids)),...
144 df_info.ndofs,df_info.ndofs);
145 reg_mat = reg_mat + A_dirichlet;
147 df_info.regularized_h10_inner_product_matrix = reg_mat;
function global_dof_index = fem_global_dof_index(params,triagrid grid)
function computing the local-to-global dof map of a fem discrete function
structure representing the fem-space information shared by all fem-functions. Implemented as handle c...