2 % structure representing the fem-space information shared by all
3 % fem-functions. Implemneted as handle class, in order to be linked
6 % B. Haasdonk 13.1.2011
8 % adapted to new assembly
16 % dimension of ranges space
19 % index vector of global DOFs
25 % number of DOFs per grid element
28 %
object of type .gridbase
31 % number of lagrange node indices
34 % the lagrange points on nodes
37 % the lagrange points on edges
40 % determinant of reference mapping
43 % global IDs of dirichlet DOFs
46 % element IDs of Robin DOFs
49 % edge IDs of Robin DOFs
52 % element IDs of Neumann DOFs
55 % edge IDs of Neumann DOFs
58 % weight matrix
for `L^2` inner product
59 l2_inner_product_matrix;
61 % weight matrix
for `H^1_0` inner product
62 h10_inner_product_matrix;
64 % weight matrix
for regularized `H^1_0` inner product
65 regularized_h10_inner_product_matrix;
70 function df_info =
feminfo (model,grid)
72 % required fields of model:
73 % pdeg: polynomial degree
74 df_info.
pdeg = model.pdeg;
75 df_info.dimrange = model.dimrange;
76 df_info.ndofs = fem_ndofs(model,grid);
77 df_info.ndofs_per_element = fem_ndofs_per_element(model);
81 df_info.lagrange_nodes = lagrange_nodes_lcoord(model.pdeg);
82 df_info.nlagrange_nodes = size(df_info.lagrange_nodes,1);
83 df_info.lagrange_edges = lagrange_nodes_edges(model.pdeg);
84 df_info.detDF = grid.A(:)*2; % determinant of reference mapping per element;
85 df_info.dirichlet_gids = [];
87 % l2_inner_product and h10_inner product
for now only
for
88 % scalar functions. Extend later.
89 % make sure, that the folliwing is called with empty
90 % dirichlet nodes, as fem_matrix_volume_component_assembly does
91 % dirichlet treatment!!!
92 % reaction term == 1 then reaction matrix is l2matrix
93 if df_info.dimrange == 1
94 model.reaction = @(grid,eindices,loc,params) ones(size(eindices,1),1);
95 model.qdeg = 2*df_info.pdeg;
96 df_info.l2_inner_product_matrix = ...
97 fem_matrix_volume_part_assembly(@fem_matrix_reac_integral_kernel,model,df_info);
98 % H10_inner_product_matrix: identity diffusion
99 model.diffusivity_tensor = @(grid,eindices,loc,params) ...
100 [ones(size(eindices,1),1), zeros(size(eindices,1),1), ...
101 zeros(size(eindices,1),1), ones(size(eindices,1),1)];
102 df_info.h10_inner_product_matrix = ...
103 fem_matrix_volume_part_assembly(@fem_matrix_diff_integral_kernel,model,df_info);
106 % find dirichlet indices as preparation
107 is_dirichlet_gid = zeros(df_info.ndofs,1);
108 % dirichlet-treatment: insert unit-vector row into A in dirichlet
109 % rows. Little loop over local lagrange node index
110 % this is repeated in r assembly as these will be kept separated
111 for i=1:df_info.nlagrange_nodes
112 % search all elements, in which the i-th Lagrange node is
114 eids = find(df_info.lagrange_edges(i,:)==1);
116 for e = eids; % possibly one point belongs to 2 edges!!
117 elids = find(grid.NBI(:,e)==-1);
118 gids = df_info.global_dof_index(elids,i);
119 is_dirichlet_gid(gids) = 1;
123 df_info.dirichlet_gids = find(is_dirichlet_gid);
125 % find robin boundary elements and edge numbers
126 ind = find(grid.NBI==-3);
127 [df_info.robin_elind,df_info.robin_edgeind] = ind2sub(size(grid.NBI),ind);
129 % find neumann boundary elements and edge numbers
130 ind = find(grid.NBI==-2);
131 [df_info.neumann_elind,df_info.neumann_edgeind] = ind2sub(size(grid.NBI),ind);
132 % regularized_h10_inner_product_matrix: make unit rows in
135 if df_info.dimrange == 1
136 reg_mat = df_info.h10_inner_product_matrix;
137 if ~isempty(df_info.dirichlet_gids)
138 reg_mat(df_info.dirichlet_gids,:) = 0;
139 reg_mat(:,df_info.dirichlet_gids) = 0;
140 A_dirichlet = sparse(df_info.dirichlet_gids,...
141 df_info.dirichlet_gids, ...
142 ones(size(df_info.dirichlet_gids)),...
143 df_info.ndofs,df_info.ndofs);
144 reg_mat = reg_mat + A_dirichlet;
146 df_info.regularized_h10_inner_product_matrix = reg_mat;