1\def\vL{\vec{L}} 2 3\section{Tools for the assemblage of linear systems}% 4\label{S:ass_tools}% 5\idx{assemblage tools|(} 6 7This section describes data structures and subroutines for matrix and 8vector assembly. Section \ref{S:update} presents basic routines for 9the update of global matrices and vectors by adding contributions from 10one single element. Data structures and routines for global matrix 11assembly are described in Section \ref{S:matrix_assemblage}. This 12includes library routines for the efficient implementation of a 13general second order linear elliptic operator. Section 14\ref{S:ass_info} presents data structures and routines for the 15handling of pre--computed integrals, which are used to speed up 16calculations in the case of problems with constant coefficients. The 17assembly of (right hand side) vectors is described in Section 18\ref{S:vector_update}. The incorporation of Dirichlet boundary values into 19the right hand side is presented in Section \ref{S:dirichlet_bound}. 20Finally, routines for generation of interpolation coefficients are 21described in Section \ref{S:I_FES}. 22 23\subsection{Element matrices and vectors}% 24\label{S:update} 25 26The usual way to assemble the system matrix and the load vector is to 27loop over all (leaf) elements, calculate the local element 28contributions and add these to the global system matrix and the global 29load vector. The updating of the load vector is rather easy. The 30contribution of a local degree of freedom is added to the value of the 31corresponding global degree of freedom. Here we have to use the 32function $j_S$ defined on each element $S$ in \mathref{book:E:global-lokal} 33on page \pageref{book:E:global-lokal}. It combines uniquely the 34local DOFs with the global ones. The basis functions provide in the 35\code{BAS\_FCTS} structure the entry \code{get\_dof\_indices()} 36which is an implementation of $j_S$, see \secref{S:basfct_data}. 37 38The updating of the system matrix is not that easy. As mentioned 39in \secref{book:S:FE-dis-2nd}, the system matrix is usually sparse and we use 40special data structures for storing these matrices, compare 41\secref{S:DOF_MATRIX}. For sparse matrices we do not have 42for each DOF a matrix row storing values for all other DOFs; only the 43values for pairs of DOFs are stored, where the corresponding 44\emph{global} basis functions have a common support. Usually, the 45exact number of entries in one row of a sparse matrix is not know a 46priori and can change during grid modifications. 47 48Thus, we use the following concept: A call of 49\code{clear\_dof\_matrix()} will not set all matrix entries to zero, 50but will remove all matrix rows from the matrix, compare the 51description of this function on page \pageref{Func:clear_dof_matrix}. 52During the updating of a matrix for the value corresponding to a pair 53of local DOFs $(i,j)$, we look in the $j_S(i)$th row of the 54matrix for a column $j_S(j)$ (the \code{col} member of 55\code{matrix\_row}); if such an entry exists, we add 56the current contribution; if this entry does not yet exist we will 57create a new entry, set the current value and column number. This 58creation may include an enlargement of the row, by linking a new 59matrix row to the list of matrix rows, if no space for a new entry is 60left. After the assemblage we then have a sparse matrix, storing all 61values for pairs of global basis functions with common support. 62 63The functions which we describe now allows also to handle matrices 64where the DOFs indexing the rows can differ from the DOFs indexing the 65columns; this makes the combination of DOFs from different finite 66element spaces possible. 67 68\begin{compatibility} 69 \label{compat:DOWB_MATRIX} 70 Previous versions of \ALBERTA defined extra-types for vector-valued 71 problems, like \code{DOF\_DOWB\_MATRIX}, \code{DOWB\_OPERATOR\_INFO} 72 etc. The ``\code{DOWB}'' (``DimOfWorldBlocks'') variants, however, 73 already incorporated all the functionality of the ordinary 74 scalar-only versions. Therefore the scalar-ony versions of most 75 data-structures have been abandoned and were replaced by the 76 ``\code{DOWB}'' variants, which in turn were renamed to use the 77 scalar-only names. For example, in the current implementation a 78 \code{DOF\_MATRIX} is in fact what older versions called a 79 \code{DOF\_DOWB\_MATRIX}; and implements the scalar-only case as 80 well as the block-matrix case. 81\end{compatibility} 82 83\subsubsection{Element matrix and vector structures} 84\label{S:elvecmat} 85 86\paragraph{Block-matrix types} 87\idx{assemblage tools!MATENT_TYPE@{\code{MATENT\_TYPE}}} 88\ddx{MATENT_TYPE@{\code{MATENT\_TYPE}}} 89\bv\begin{lstlisting}[label=T:MATENT_TYPE,name=MATENT_TYPE,caption={[data-type: \code{MATENT\_TYPE}]}] 90typedef enum matent_type { 91 MATENT_NONE = -1, 92 MATENT_REAL = 0, 93 MATENT_REAL_D = 1, 94 MATENT_REAL_DD = 2 95} MATENT_TYPE; 96\end{lstlisting}\ev 97 98\noindent 99Description: This enumeration type defines symbolic types for 100block-matrix entries. \code{MATENT\_REAL} means scalar blocks, 101\code{MATENT\_REAL\_D} stands for diagonal blocks, and 102\code{MATENT\_REAL\_DD} is a code for full matrix blocks. In general, 103data-structures make use of these types to store the matrix blocks in 104an efficient way. 105 106\paragraph{Structure for element matrices} 107\idx{assemblage tools!EL_MATRIX@{\code{EL\_MATRIX}}} 108\ddx{EL_MATRIX@{\code{EL\_MATRIX}}} 109\bv\begin{lstlisting}[label=T:EL_MATRIX,name=EL_MATRIX,caption={[data-type: \code{EL\_MATRIX}]}] 110typedef struct el_matrix EL_MATRIX; 111 112struct el_matrix 113{ 114 MATENT_TYPE type; 115 int n_row, n_col; 116 int n_row_max, n_col_max; 117 union { 118 REAL *const*real; 119 REAL_D *const*real_d; 120 REAL_DD *const*real_dd; 121 } data; 122 DBL_LIST_NODE row_chain; 123 DBL_LIST_NODE col_chain; 124}; 125\end{lstlisting}\ev 126 127\noindent 128Description: A data structure to store per-element contributions 129during the assembling of discrete systems. There is some limited 130support for the operation of element-matrices on element-vectors and 131global DOF-vectors, see \secref{S:elementblas}. 132\begin{descr} 133 \kitem{type} One out of \code{MATENT\_REAL}, \code{MATENT\_REAL\_D} 134 or \code{MATENT\_REAL\_DD}. The entries stored in 135 \code{EL\_MATRIX->data} have to be interpreted accordingly. See 136 \code{MATENT\_TYPE} on page \pageref{T:MATENT_TYPE}. 137%% 138 \kitem{n\_row} is the number of rows of the element matrix 139%% 140 \kitem{n\_col} is the number of columns of the element matrix 141%% 142 \kitem{n\_row\_max} is the maximal number of rows. The number of rows 143 can vary from element to element if the underlying basis functions 144 have a per-element initializer. 145%% 146 \kitem{n\_col\_max} is the maximal number of columns. 147% 148 \kitem{data, data.real, data.real\_d, data.real\_dd} 149 \code{EL\_MATRIX->data} is a union, its components should be accessed 150 according to the symmetry type indicated by \code{EL\_MATRIX->type}. 151 %% 152 \kitem{row\_chain, col\_chain} If the underlying finite element 153 spaces are a direct sum of function spaces, then the resulting 154 element matrices have a block-layout. The link to the other parts of 155 the resulting block-matrix is implemented using cyclic doubly linked 156 lists, \code{row\_chain} and \code{col\_chain} are the corresponding 157 list-nodes. There is a separate section explaining how to handle such 158 chains of objects, see \secref{S:chain_impl}. 159\end{descr} 160 161\paragraph{Structures for element vectors} 162\ddx{EL_INT_VEC@{\code{EL\_INT\_VEC}}} 163\idx{assemblage tools!EL_INT_VEC@{\code{EL\_INT\_VEC}}} 164\ddx{EL_DOF_VEC@{\code{EL\_DOF\_VEC}}} 165\idx{assemblage tools!EL_DOF_VEC@{\code{EL\_DOF\_VEC}}} 166\ddx{EL_UCHAR_VEC@{\code{EL\_UCHAR\_VEC}}} 167\idx{assemblage tools!EL_UCHAR_VEC@{\code{EL\_UCHAR\_VEC}}} 168\ddx{EL_SCHAR_VEC@{\code{EL\_SCHAR\_VEC}}} 169\idx{assemblage tools!EL_SCHAR_VEC@{\code{EL\_SCHAR\_VEC}}} 170\ddx{EL_BNDRY_VEC@{\code{EL\_BNDRY\_VEC}}} 171\idx{assemblage tools!EL_BNDRY_VEC@{\code{EL\_BNDRY\_VEC}}} 172\ddx{EL_PTR_VEC@{\code{EL\_PTR\_VEC}}} 173\idx{assemblage tools!EL_PTR_VEC@{\code{EL\_PTR\_VEC}}} 174\ddx{EL_REAL_VEC@{\code{EL\_REAL\_VEC}}} 175\idx{assemblage tools!EL_REAL_VEC@{\code{EL\_REAL\_VEC}}} 176\ddx{EL_REAL_VEC_D@{\code{EL\_REAL\_VEC\_D}}} 177\idx{assemblage tools!EL_REAL_VEC_D@{\code{EL\_REAL\_VEC\_D}}} 178\ddx{EL_REAL_D_VEC@{\code{EL\_REAL\_D\_VEC}}} 179\idx{assemblage tools!EL_REAL_D_VEC@{\code{EL\_REAL\_D\_VEC}}} 180\bv\begin{lstlisting}[label=T:EL_INT_VEC,name=EL_XXX_VEC,caption={[data-types: \code{EL\_XXX\_VEC}]}] 181typedef struct el_int_vec EL_INT_VEC; 182typedef struct el_dof_vec EL_DOF_VEC; 183typedef struct el_uchar_vec EL_UCHAR_VEC; 184typedef struct el_schar_vec EL_SCHAR_VEC; 185typedef struct el_bndry_vec EL_BNDRY_VEC; 186typedef struct el_ptr_vec EL_PTR_VEC; 187typedef struct el_real_vec EL_REAL_VEC; 188typedef struct el_real_vec_d EL_REAL_VEC_D; 189typedef struct el_real_d_vec EL_REAL_D_VEC; 190\end{lstlisting}\ev 191\label{T:EL_DOF_VEC} 192\label{T:EL_UCHAR_VEC} 193\label{T:EL_SCHAR_VEC} 194\label{T:EL_BNDRY_VEC} 195\label{T:EL_PTR_VEC} 196\label{T:EL_REAL_D_VEC} 197The \code{el\_*\_vec} structures are declared similary, the only 198difference between them is the type of the structure entry vec. Below, 199the \code{EL\_REAL\_VEC} structure is given: 200%% 201\bv\begin{lstlisting}[label=T:EL_REAL_VEC,name=EL_REAL_VEC,caption={data-type: \code{EL\_REAL\_VEC}}] 202struct el_real_vec 203{ 204 int n_components; 205 int n_components_max; 206 DBL_LIST_NODE chain; 207 int reserved; 208 REAL vec[1]; /* different type in EL_INT_VEC, ... */ 209}; 210\end{lstlisting}\ev 211%% 212and the \code{EL\_REAL\_VEC\_D} structure is described in detail: 213%% 214\bv\begin{lstlisting}[label=T:EL_REAL_VEC_D,name=EL_REAL_VEC_D] 215struct el_real_vec_d 216{ 217 int n_components; 218 int n_components_max; 219 DBL_LIST_NODE chain; 220 int stride; /* either 1 or DIM_OF_WORLD */ 221 REAL vec[1]; 222}; 223\end{lstlisting}\ev 224%% 225Description: 226\begin{descr} 227 \kitem{n\_components} The actual number of components available in 228 and following \code{EL\_XXX\_VEC->vec}. Note that the actual number 229 of components is -- of course -- larger than $1$ in general 230 \code{get\_el\_XXX\_vec(bas\_fcts)} takes care of allocating enough 231 space. 232% 233 \kitem{n\_components\_max} Behing \code{EL\_XXX\_VEC->vec[0]} may 234 actually be more space available than the number of currently valid 235 entries as indicated by \code{EL\_XXX\_VEC->n\_components}; this is 236 the maximum size to access without crossing the bounds of the data 237 segment allocated for this element vector. 238% 239 \kitem{chain} If the underlying basis-function implementation is 240 part of a chain of sets of basis functions, then this structure is 241 inherited also by the element vectors: they are chained using a 242 doubly linked list, \code{chain} is the corresponding list-node. 243 There is a separate section about such chained objects, see 244 \secref{S:chain_impl}. 245% 246 \kitem{stride, reserved} For element vectors other than an 247 \code{EL\_REAL\_VEC\_D} this is a reserved value and actually tied to 248 the constant value $1$ with the exception of a 249 \code{EL\_REAL\_D\_VEC} were \code{reserved} is fixed at 250 \code{DIM\_OF\_WORLD}. For \code{EL\_REAL\_VEC\_D} this varies, 251 based on the dimension of the range of the underlying basis function 252 implementation. For vector-valued basis functions 253 \code{EL\_REAL\_VEC\_D->stride} is again tied to $1$, for 254 scalar-valued basis functions \code{EL\_REAL\_VEC\_D->stride} is 255 fixed at \code{DIM\_OF\_WORLD}, in both cases it gives the number of 256 \code{REAL}'s belonging to a single \code{DOF}. See also 257 \code{DOF\_REAL\_VEC\_D} on page \pageref{T:DOF_REAL_VEC_D}. 258% 259 \kitem{vec[1]} Start of the data-segment, 260 \code{EL\_XXX\_VEC->n\_components} items contain valid data, 261 \code{EL\_XXX\_VEC->n\_components\_max} items are allocated. Note 262 that for a \code{EL\_REAL\_VEC\_D} vector the numbers have to be 263 multiplied by \code{EL\_REAL\_VEC\_D->stride} to get the actual 264 number of \code{REAL}'s allocated. 265\end{descr} 266 267\subsubsection{Accumulating per-element contributions} 268The following functions can be used on elements for updating matrices 269and vectors. 270\fdx{add_element_matrix()@{\code{add\_element\_matrix()}}} 271\idx{assemblage tools!add_element_matrix()@{\code{add\_element\_matrix()}}} 272\fdx{add_element_vec()@{\code{add\_element\_vec()}}} 273\idx{assemblage tools!add_element_vec()@{\code{add\_element\_vec()}}} 274\fdx{add_element_d_vec()@{\code{add\_element\_d\_vec()}}} 275\idx{assemblage tools!add_element_d_vec()@{\code{add\_element\_d\_vec()}}} 276\fdx{add_element_vec_dow()@{\code{add\_element\_vec\_dow()}}} 277\idx{assemblage tools!add_element_vec_dow()@{\code{add\_element\_vec\_dow()}}} 278\bv\begin{lstlisting}[name={proto-type add_element_matrix()},label=C:add_element_matrix] 279void add_element_matrix(DOF_MATRIX *matrix, REAL factor, 280 const EL_MATRIX *el_matrix, MatrixTranspose transpose, 281 const EL_DOF_VEC *row_dof, const EL_DOF_VEC *col_dof, 282 const EL_SCHAR_VEC *bound); 283void add_element_vec(DOF_REAL_VEC *drv, REAL factor, const EL_REAL_VEC *el_vec, 284 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 285void add_element_d_vec(DOF_REAL_D_VEC *drdv, REAL factor, 286 const EL_REAL_D_VEC *el_vec, const EL_DOF_VEC *dof, 287 const EL_SCHAR_VEC *bound); 288void add_element_vec_dow(DOF_REAL_VEC_D *drdv, REAL factor, 289 const EL_REAL_VEC_D *el_vec, const EL_DOF_VEC *dof, 290 const EL_SCHAR_VEC *bound); 291\end{lstlisting}\ev 292\paragraph{Descriptions} 293\begin{descr} 294\kitem{add\_element\_matrix(mat, 295\fdx{add_element_matrix()@{\code{add\_element\_matrix()}}} 296\idx{assemblage tools!add_element_matrix()@{\code{add\_element\_matrix()}}} 297 factor, el\_mat, transpose, row\_dof, col\_dof, bound)}\hfill 298 299Updates the global \code{DOF\_MATRIX} \code{mat} by adding element 300contributions. If \code{row\_dof} equals \code{col\_dof}, the diagonal 301element is \emph{always} the first entry in a matrix row; this makes 302the access to the diagonal element easy for a diagonal preconditioner, 303for example. In general, \code{add\_element\_matrix()} does the 304following: for all \code{i} the values 305\code{fac*el\_mat->data.\{REAL,REAL\_D,REAL\_DD\}[i][j]} are added to 306the entries at the position (\code{row\_dof->vec[i],col\_dof->vec[j]}) 307in the global matrix \code{mat} 308($\code{0}\leq\code{i}<\code{el\_mat->n\_row}$, 309$\code{0}\leq\code{j}<\code{el\_mat->n\_col}$). If such an entry 310exists in the row number \code{row\_dof->vec[i]} the global matrix 311\code{mat} the value is simply added. Otherwise a new entry is created 312in the row, the value is set and the column number is set to 313\code{col\_dof[j]}. This may imply an enlargement of the row by adding 314a new \code{MATRIX\_ROW} structure to the list of matrix rows. 315 316Note that the first element matrix added to \code{mat} after calling 317\code{clear\_dof\_matrix()} determines the block-type of the global 318matrix \code{mat}. It is possible to add element-matrices with higher 319block-symmetry to global \code{DOF\_MATRIX}es with lower 320block-symmetry, for example it is allowed to add \code{el\_mat} to 321\code{mat} if \code{el\_mat->type == MATENT\_REAL} and \code{mat->type 322 == MATENT\_REAL\_DD}. 323 324\begin{description} 325\item[Parameters]\hfill 326 \begin{description} 327 \item[\code{mat}] the global \code{DOF\_MATRIX}. 328 \item[\code{factor}] is a multiplier for the element contributions; 329 usually \code{factor} is \code{1} or \code{-1}; 330 331 \item[\code{el\_mat}] is a matrix of size 332 $\code{n\_row}\times\code{n\_col}$ storing the element 333 contributions; 334 335 \item[\code{transpose}] the original matrix is used if 336 \code{transpose} == \code{NoTranspose} (= 0) and the transposed 337 matrix if \code{transpose} == \code{Transpose} (= 1); 338 339 \item[\code{row\_dof}] is a vector of length 340 \code{row\_dof->n\_components} storing the global row indices; 341 342 \item[\code{col\_dof}] is a vector of length 343 \code{col\_dof->n\_components} storing the global column indices, 344 \code{col\_dof} may be a \nil pointer if the DOFs indexing the 345 columns are the same as the DOFs indexing the rows; in this case 346 \code{col\_dof = row\_dof} is used; 347 348 \item[\code{bound}] is either \code{NULL} or an 349 \code{EL\_SCHAR\_Vec} stucture storing a vector of length 350 \code{bound->n\_components}. In this case 351 \code{bound->n\_components} must match either 352 \code{row\_dof->n\_components} or \code{col\_dof->n\_components}, 353 depending on the value of \code{transpose}. 354 355 If \code{bound->vec[i] >= DIRICHLET}, then the following happens: 356 \begin{description} 357 \item[\code{row\_dof == col\_dof}] In the global \code{matrix} the 358 row \code{row\_dof->vec[i]} is cleared to zero, with the 359 exception of the diagonal entry, which is set to \code{1.0}. 360 \item[\code{row\_dof != col\_dof}] In the global \code{matrix} the 361 row \code{row\_dof->vec[i]} is cleared to zero. 362 \end{description} 363 364 All other contributions of \code{el\_mat} are added to 365 \code{matrix} as usual. This allows for a convenient way to 366 implement inhomogeneous Dirichlet boundary conditions, without 367 having to modify the right-hand-side of the discrete systems 368 explicitly. 369 \end{description} 370\end{description} 371 372 \hrulefill 373\fdx{add_element_vec()@{\code{add\_element\_vec()}}} 374\idx{assemblage tools!add_element_vec()@{\code{add\_element\_vec()}}} 375\fdx{add_element_d_vec()@{\code{add\_element\_d\_vec()}}} 376\idx{assemblage tools!add_element_d_vec()@{\code{add\_element\_d\_vec()}}} 377\fdx{add_element_vec_dow()@{\code{add\_element\_vec\_dow()}}} 378\idx{assemblage tools!add_element_vec_dow()@{\code{add\_element\_vec\_dow()}}} 379 \kitem{add\_element\_vec(drv, factor, el\_vec, dof, bound)} 380 \kitem{add\_element\_d\_vec(drv, factor, el\_vec, dof, bound)} 381 \kitem{add\_element\_vec\_d(drv, factor, el\_vec, dof, bound)}\hfill 382 383 These do similar things as \code{add\_element\_matrix()}, but with 384 element vectors. \secref{S:elementblas} also lists other routines 385 which might be helpful in this context. 386\end{descr} 387 388\subsubsection{Allocation and filling of element vectors} 389\label{S:fillgetelvec} 390 391\paragraph{Prototypes} 392\fdx{get_dof_indices()@{\code{get\_dof\_indices()}}} 393\fdx{get_bound()@{\code{get\_bound()}}} 394\fdx{el_interpol()@{\code{el\_interpol()}}} 395\fdx{el_interpol_dow()@{\code{el\_interpol\_dow()}}} 396\fdx{dirichlet_map()@{\code{dirichlet\_map()}}} 397%% 398%\fdx{init_el_int_vec()@{\code{init\_el\_int\_vec()}}} 399%\fdx{init_el_dof_vec()@{\code{init\_el\_dof\_vec()}}} 400%\fdx{init_el_real_vec()@{\code{init\_el\_real\_vec()}}} 401%\fdx{init_el_real_d_vec()@{\code{init\_el\_real\_d\_vec()}}} 402%\fdx{init_el_uchar_vec()@{\code{init\_el\_uchar\_vec()}}} 403%\fdx{init_el_schar_vec()@{\code{init\_el\_schar\_vec()}}} 404%\fdx{init_el_ptr_vec()@{\code{init\_el\_ptr\_vec()}}} 405%% 406\fdx{fill_el_int_vec()@{\code{fill\_el\_int\_vec()}}} 407\fdx{fill_el_real_vec()@{\code{fill\_el\_real\_vec()}}} 408\fdx{fill_el_real_d_vec()@{\code{fill\_el\_real\_d\_vec()}}} 409\fdx{fill_el_real_vec_d()@{\code{fill\_el\_real\_vec\_d()}}} 410\fdx{fill_el_uchar_vec()@{\code{fill\_el\_uchar\_vec()}}} 411\fdx{fill_el_schar_vec()@{\code{fill\_el\_schar\_vec()}}} 412%% 413\fdx{get_el_int_vec()@{\code{get\_el\_int\_vec()}}} 414\fdx{get_el_dof_vec()@{\code{get\_el\_dof\_vec()}}} 415\fdx{get_el_uchar_vec()@{\code{get\_el\_uchar\_vec()}}} 416\fdx{get_el_schar_vec()@{\code{get\_el\_schar\_vec()}}} 417\fdx{get_el_bndry_vec()@{\code{get\_el\_bndry\_vec()}}} 418\fdx{get_el_ptr_vec()@{\code{get\_el\_ptr\_vec()}}} 419\fdx{get_el_real_vec()@{\code{get\_el\_real\_vec()}}} 420\fdx{get_el_real_d_vec()@{\code{get\_el\_real\_d\_vec()}}} 421\fdx{get_el_real_vec_d()@{\code{get\_el\_real\_vec\_d()}}} 422%% 423\fdx{free_el_int_vec()@{\code{free\_el\_int\_vec()}}} 424\fdx{free_el_dof_vec()@{\code{free\_el\_dof\_vec()}}} 425\fdx{free_el_uchar_vec()@{\code{free\_el\_uchar\_vec()}}} 426\fdx{free_el_schar_vec()@{\code{free\_el\_schar\_vec()}}} 427\fdx{free_el_bndry_vec()@{\code{free\_el\_bndry\_vec()}}} 428\fdx{free_el_ptr_vec()@{\code{free\_el\_ptr\_vec()}}} 429\fdx{free_el_real_vec()@{\code{free\_el\_real\_vec()}}} 430\fdx{free_el_real_d_vec()@{\code{free\_el\_real\_d\_vec()}}} 431\fdx{free_el_real_vec_d()@{\code{free\_el\_real\_vec\_d()}}} 432\bv\begin{lstlisting}[label=P:ELVEC_PROTOS,name=ELVEC_PROTOS,caption={[proto-types: filling element vectors]}] 433EL_DOF_VEC *get_dof_indices(EL_DOF_VEC *dofs, const FE_SPACE *fe_space, 434 const EL *el); 435EL_BNDRY_VEC *get_bound(EL_BNDRY_VEC *bndry, const BAS_FCTS *bas_fcts, 436 const EL_INFO *el_info); 437void el_interpol(EL_REAL_VEC *coeff, const EL_INFO *el_info, int wall, 438 const EL_INT_VEC *indices, LOC_FCT_AT_QP f, void *ud, 439 const BAS_FCTS *bas_fcts); 440void el_interpol_dow(EL_REAL_VEC_D *coeff, const EL_INFO *el_info, int wall, 441 const EL_INT_VEC *indices, LOC_FCT_D_AT_QP f, 442 void *f_data, const BAS_FCTS *bas_fcts); 443void dirichlet_map(EL_SCHAR_VEC *bound, const EL_BNDRY_VEC *bndry_bits, 444 const BNDRY_FLAGS mask); 445 446const EL_INT_VEC * 447fill_el_int_vec(EL_INT_VEC *el_vec, EL *el, const DOF_INT_VEC *dof_vec); 448const EL_REAL_VEC * 449fill_el_real_vec(EL_REAL_VEC *el_vec, EL *el, const DOF_REAL_VEC *dof_vec); 450const EL_REAL_D_VEC * 451fill_el_real_d_vec(EL_REAL_D_VEC *el_vec, EL *el, const DOF_REAL_D_VEC *dof_vec); 452const EL_REAL_VEC_D * 453fill_el_real_vec_d(EL_REAL_VEC_D *el_vec, EL *el, const DOF_REAL_VEC_D *dof_vec); 454const EL_UCHAR_VEC * 455fill_el_uchar_vec(EL_UCHAR_VEC *el_vec, EL *el, const DOF_UCHAR_VEC *dof_vec); 456const EL_SCHAR_VEC * 457fill_el_schar_vec(EL_SCHAR_VEC *el_vec, EL *el, const DOF_SCHAR_VEC *dof_vec); 458 459EL_INT_VEC *get_el_int_vec(const BAS_FCTS *bas_fcts); 460EL_DOF_VEC *get_el_dof_vec(const BAS_FCTS *bas_fcts); 461EL_UCHAR_VEC *get_el_uchar_vec(const BAS_FCTS *bas_fcts); 462EL_SCHAR_VEC *get_el_schar_vec(const BAS_FCTS *bas_fcts); 463EL_BNDRY_VEC *get_el_bndry_vec(const BAS_FCTS *bas_fcts); 464EL_PTR_VEC *get_el_ptr_vec(const BAS_FCTS *bas_fcts); 465EL_REAL_VEC *get_el_real_vec(const BAS_FCTS *bas_fcts); 466EL_REAL_D_VEC *get_el_real_d_vec(const BAS_FCTS *bas_fcts); 467EL_REAL_VEC_D *get_el_real_vec_d(const BAS_FCTS *bas_fcts); 468 469void free_el_int_vec(EL_INT_VEC *el_vec); 470void free_el_dof_vec(EL_DOF_VEC *el_vec); 471void free_el_uchar_vec(EL_UCHAR_VEC *el_vec); 472void free_el_schar_vec(EL_SCHAR_VEC *el_vec); 473void free_el_bndry_vec(EL_BNDRY_VEC *el_vec); 474void free_el_ptr_vec(EL_PTR_VEC *el_vec); 475void free_el_real_vec(EL_REAL_VEC *el_vec); 476void free_el_real_d_vec(EL_REAL_D_VEC *el_vec); 477void free_el_real_vec_d(EL_REAL_VEC_D *el_vec); 478 479DEF_EL_VEC_VAR(VECNAME, name, _size, _size_max, _init); 480DEF_EL_VEC_CONST(VECNAME, name, _size, _size_max); 481ALLOC_EL_VEC(VECNAME, _size, _size_max); 482\end{lstlisting}\ev 483 484\paragraph{Descriptions} 485\begin{descr} 486 \kitem{get\_dof\_indices(dofs, fe\_space, el)} 487 488 Compute the mapping 489 \fdx{get_dof_indices()@{\code{get\_dof\_indices()}}} 490 between the local \code{DOF}-indices on \code{el} and the global 491 \code{DOF}-indices according to \code{fe\_space->admin}. 492 \begin{description} 493 \item[Parameters]~ 494 \begin{description} 495 \item[\code{dofs}] Storage for the result or \code{NULL}. In the 496 latter case the mapping is returned in a statically allocated 497 \code{EL\_DOF\_VEC}. \emph{Note: this storage area will be 498 overwritten on the next call to this function, even if the 499 \code{fe\_space} argument differs.} 500 \item[\code{fe\_space}] The finite element space to compute the mapping for. 501 \item[\code{el}] The current mesh element (\emph{not} the current 502 \code{EL\_INFO} pointer, use \code{EL\_INFO->el}). 503 \end{description} 504 \item[return] Either again the argument \code{dofs} or -- if 505 \code{dofs == NULL} -- a pointer to a statically allocated 506 \code{EL\_DOF\_VEC}. 507 \item[examples] 508\begin{samepage} 509With pre-allocated \code{EL\_DOF\_VEC}: 510\bv\begin{lstlisting}[caption={[Example for \code{get\_dof\_indices()}]},label=C:get_dof_indices_example] 511EL_DOF_VEC *dofs = get_el_dof_vec(fe_space->bas_fcts); 512TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL) { 513 int i; 514 515 get_dof_indices(dofs, fe_space, el_info->el); 516 for (i = 0; i < bas_fcts->n_bas_fcts; i++) { 517 MSG("dofs[%d] = %d\n", dofs->vec[i]); 518 } 519} TRAVERSE_NEXT(); 520free_el_dof_vec(dofs); 521\end{lstlisting}\ev 522\end{samepage} 523 524\begin{samepage} 525Without pre-allocated \code{EL\_DOF\_VEC}: 526\bv\begin{lstlisting}[caption={[Example for \code{get\_dof\_indices()}]},label=C:get_dof_indices_example2] 527TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL) { 528 int i; 529 EL_DOF_VEC *dofs = get_dof_indices(NULL, fe_space, el_info->el); 530 531 for (i = 0; i < bas_fcts->n_bas_fcts; i++) { 532 MSG("dofs[%d] = %d\n", dofs->vec[i]); 533 } 534} TRAVERSE_NEXT(); 535\end{lstlisting}\ev 536\end{samepage} 537%%\item[see also] 538 \end{description} 539 %% 540 \hrulefill 541 %% 542 \kitem{get\_bound(bndry, bas\_fcts, el\_info)} 543 \fdx{get_bound()@{\code{get\_bound()}}} Extract the boundary types of 544 the local DOFs of \code{bas\_fcts}. The boundary types are returned 545 in form of a bit-mask. If bit \code{j} in the bit-mask 546 \code{bndry[i]} is set, then the local DOF number \code{i} belongs to 547 the boundary segment which has been assigned the number \code{j} in 548 the macro-triangulation. Boundary types range from $1$ to $255$. 549 \begin{description} 550 \item[Parameters]~ 551 \begin{description} 552 \item[\code{EL\_BNDRY\_VEC *bndry}] Storage for the result or 553 \code{NULL}. In the latter case the data is returned in a 554 statically allocated \code{EL\_BNDRY\_VEC}. 555 \item[\code{BAS\_FCTS *bas\_fcts}] The local basis functions. 556 \item[\code{const EL\_INFO *el\_info}] The current mesh element 557 info structure. (\emph{not} the current \code{EL\_INFO} pointer. 558 \end{description} 559 \item[return] Either again the argument \code{bndry} or -- if 560 \code{bndry == NULL} -- a pointer to a statically allocated 561 \code{EL\_BNDRY\_VEC}. 562 \item[examples] 563\begin{samepage} 564With pre-allocated \code{EL\_BNDRY\_VEC}: 565 \bv\begin{lstlisting} 566EL_BNDRY_VEC *bndry = get_el_bndry_vec(bas_fcts); 567TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_BOUND) { 568 int i, j; 569 get_bound(bndry, bas_fcts, el_info); 570 for (i = 0; i < bas_fcts->n_bas_fcts; i++) { 571 for (j = 1; j < N_BNDRY_TYPES; j++) { 572 if (BNNDRY_FLAGS_IS_INTERIOR(bndry->vec[i])) { 573 MSG("Local dof %d is an interior DOF\n"); 574 } else if (BNDRY_FLAGS_IS_AT_BNDRY(bndry->vec[i], j)) { 575 MSG("Local dof %d belongs to boundary segment %d\n", i, j); 576 } 577 } 578 } 579} TRAVERSE_NEXT(); 580free_el_bndry_vec(bndry); 581 \end{lstlisting}\ev 582\end{samepage} 583 584\begin{samepage} 585Without pre-allocated \code{EL\_BNDRY\_VEC}: 586 \bv\begin{lstlisting} 587TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL) { 588 int i,j; 589 EL_BNDRY_VEC *bndry = get_bound(NULL, bas_fcts, el_info); 590 591 for (i = 0; i < bas_fcts->n_bas_fcts; i++) { 592 for (j = 1; j < N_BNDRY_TYPES; j++) { 593 if (BNNDRY_FLAGS_IS_INTERIOR(bndry->vec[i])) { 594 MSG("Local dof %d is an interior DOF\n"); 595 } else if (BNDRY_FLAGS_IS_AT_BNDRY(bndry->vec[i], j)) { 596 MSG("Local dof %d belongs to boundary segment %d\n", i, j); 597 } 598 } 599 } 600} TRAVERSE_NEXT(); 601 \end{lstlisting}\ev 602\end{samepage} 603 \end{description} 604 %% 605 \hrulefill 606 %% 607 \kitem{fill\_el\_int\_vec(el\_vec, el, dof\_vec)} 608 \kitem{fill\_el\_real\_vec(el\_vec, el, dof\_vec)} 609 \kitem{fill\_el\_real\_d\_vec(el\_vec, el, dof\_vec)} 610 \kitem{fill\_el\_real\_vec\_d(el\_vec, el, dof\_vec)} 611 \kitem{fill\_el\_uchar\_vec(el\_vec, el, dof\_vec)} 612 \kitem{fill\_el\_schar\_vec(el\_vec, el, dof\_vec)} 613 \fdx{fill_el_int_vec()@{\code{fill\_el\_int\_vec()}}} 614 \fdx{fill_el_real_vec()@{\code{fill\_el\_real\_vec()}}} 615 \fdx{fill_el_real_d_vec()@{\code{fill\_el\_real\_d\_vec()}}} 616 \fdx{fill_el_real_vec_d()@{\code{fill\_el\_real\_vec\_d()}}} 617 \fdx{fill_el_uchar_vec()@{\code{fill\_el\_uchar\_vec()}}} 618 \fdx{fill_el_schar_vec()@{\code{fill\_el\_schar\_vec()}}}~\hfill 619 %% 620 621 Fill the respective element vector with data. The description below 622 is for \code{fill\_el\_real\_vec()}, the other versions work similar. 623 \begin{description} 624 \item[Parameters]\hfill 625 \begin{description} 626 \item[\code{EL\_REAL\_VEC *el\_vec}] Storage for the result or 627 \code{NULL}. In the latter case the return value is 628 \hyperlink{DOF_REAL_VEC:vec_loc}{\code{DOF\_REAL\_VEC->vec\_loc}}; 629 the data will be overwritten on the next call to 630 \code{fill\_el\_real\_vec()} with the same \code{dof\_vec} 631 argument. Calling \code{fill\_el\_real\_vec()} with \emph{other} 632 DOF-vectors will \emph{not} invalidate the data. 633 \item[\code{EL *el}] The current mesh element (\emph{not} the 634 current \code{EL\_INFO} pointer, use \code{EL\_INFO->el}). 635 \item[\code{DOF\_REAL\_VEC *dof\_vec}] The global \code{DOF}-vector 636 to extract the data from. 637 \end{description} 638 \item[return] Either again a the pointer \code{el\_vec} or -- 639 if \code{el\_vec == NULL} a pointer to a statically allocated 640 result space which will be overwritten on the next call to 641 \code{fill\_el\_real\_vec()}. \emph{Warning:} see ``bugs'' below. 642 \end{description} 643 %% 644 \hrulefill 645 %% 646 \kitem{get\_el\_int\_vec(bas\_fcts)} 647 \kitem{get\_el\_dof\_vec(bas\_fcts)} 648 \kitem{get\_el\_uchar\_vec(bas\_fcts)} 649 \kitem{get\_el\_schar\_vec(bas\_fcts)} 650 \kitem{get\_el\_bndry\_vec(bas\_fcts)} 651 \kitem{get\_el\_ptr\_vec(bas\_fcts)} 652 \kitem{get\_el\_real\_vec(bas\_fcts)} 653 \kitem{get\_el\_real\_d\_vec(bas\_fcts)} 654 \kitem{get\_el\_real\_vec\_d(bas\_fcts)}\hfill 655\fdx{get_el_int_vec()@{\code{get\_el\_int\_vec()}}} 656\fdx{get_el_dof_vec()@{\code{get\_el\_dof\_vec()}}} 657\fdx{get_el_uchar_vec()@{\code{get\_el\_uchar\_vec()}}} 658\fdx{get_el_schar_vec()@{\code{get\_el\_schar\_vec()}}} 659\fdx{get_el_bndry_vec()@{\code{get\_el\_bndry\_vec()}}} 660\fdx{get_el_ptr_vec()@{\code{get\_el\_ptr\_vec()}}} 661\fdx{get_el_real_vec()@{\code{get\_el\_real\_vec()}}} 662\fdx{get_el_real_d_vec()@{\code{get\_el\_real\_d\_vec()}}} 663\fdx{get_el_real_vec_d()@{\code{get\_el\_real\_vec\_d()}}} 664 665The \code{get\_el\_*\_vec()} routines automatically allocates enough 666memory for the element data vector \code{vec} as indicated by 667\code{bas\_fcts->n\_bas\_fcts}. 668\begin{description} 669\item[Parameters] 670 \begin{description} 671 \item[\code{const BAS\_FCTS *bas\_fcts}] 672 \end{description} 673\item[return] A pointer to a dynamically allocated element 674 vector of the respective type. 675\item[examples] See the first example for the 676 \code{fill\_el\_real\_vec()} function. 677\end{description} 678 %% 679 \hrulefill 680 %% 681 \kitem{free\_el\_int\_vec(el\_vec)} 682 \kitem{free\_el\_dof\_vec(el\_vec)} 683 \kitem{free\_el\_uchar\_vec(el\_vec)} 684 \kitem{free\_el\_schar\_vec(el\_vec)} 685 \kitem{free\_el\_bndry\_vec(el\_vec)} 686 \kitem{free\_el\_ptr\_vec(el\_fcts)} 687 \kitem{free\_el\_real\_vec(bas\_fcts)} 688 \kitem{free\_el\_real\_d\_vec(bas\_fcts)} 689 \kitem{free\_el\_real\_vec\_d(bas\_fcts)} 690\fdx{free_el_int_vec()@{\code{free\_el\_int\_vec()}}} 691\fdx{free_el_dof_vec()@{\code{free\_el\_dof\_vec()}}} 692\fdx{free_el_uchar_vec()@{\code{free\_el\_uchar\_vec()}}} 693\fdx{free_el_schar_vec()@{\code{free\_el\_schar\_vec()}}} 694\fdx{free_el_bndry_vec()@{\code{free\_el\_bndry\_vec()}}} 695\fdx{free_el_ptr_vec()@{\code{free\_el\_ptr\_vec()}}} 696\fdx{free_el_real_vec()@{\code{free\_el\_real\_vec()}}} 697\fdx{free_el_real_d_vec()@{\code{free\_el\_real\_d\_vec()}}} 698\fdx{free_el_real_vec_d()@{\code{free\_el\_real\_vec\_d()}}} 699\hfill 700 701The \code{free\_el\_XXX\_vec()} routines free all previously allocated 702storage for \code{el\_XXX\_vec} data. 703\begin{description} 704\item[Parameters] 705 \begin{description} 706 \item[\code{const BAS\_FCTS *bas\_fcts}] 707 \end{description} 708\item[return] \code{void} 709\item[examples] See the first example for the 710 \code{fill\_el\_real\_vec()} function. 711\end{description} 712 %% 713 \hrulefill 714 %% 715\kitem{DEF\_EL\_VEC\_VAR(VECNAME, name, size, size\_max, init)} 716\fdx{DEF_EL_VEC_VAR()@{\code{DEF\_EL\_VEC\_VAR()}}} 717\mdx{DEF_EL_VEC_VAR()@{\code{DEF\_EL\_VEC\_VAR()}}}\hfill 718 719This is a macro which defines a (local) variable with id \code{name}, 720pointing to an \code{EL\_VECNAME\_VEC} of size \code{size}, holding a 721maximal number of elements \code{max\_size}, which is initialised if 722\code{init} is \code{true}. \code{size} and \code{size\_max} may be 723variables. 724 725 \hrulefill 726 %% 727\kitem{DEF\_EL\_VEC\_CONST(VECNAME, name, size, size\_max)}\hfill 728\fdx{DEF_EL_VEC_CONST()@{\code{DEF\_EL\_VEC\_CONST()}}} 729\mdx{DEF_EL_VEC_CONST()@{\code{DEF\_EL\_VEC\_CONST()}}} 730 731This is a macro which defines a (local) variable with id \code{name}, 732pointing to an \code{EL\_VECNAME\_VEC} of size \code{size}, holding a 733maximal number of elements \code{max\_size}. \code{size} and 734\code{size\_max} must be constant values. 735 736 \hrulefill 737 %% 738\kitem{ALLOC\_EL\_VEC(VECNAME, size, size\_max)} 739\fdx{ALLOC_EL_VEC()@{\code{ALLOC\_EL\_VEC()}}} 740\mdx{ALLOC_EL_VEC()@{\code{ALLOC\_EL\_VEC()}}} 741 742This macro allocates a \code{EL\_VECNAME\_VEC} with enough storage to 743hold \code{size\_max} elements; the \code{n\_components} component of 744the element vector structure is set to \code{size}. 745 746\hrulefill 747 %% 748 \kitem{el\_interpol(coeff, el\_info, wall, indices, f, ud, bas\_fcts)} 749 \fdx{el_interpol()@{\code{el\_interpol()}}} 750 %% 751 \kitem{el\_interpol\_dow(coeff, el\_info, wall, indices, f, f\_data, ud, bas\_fcts)} 752 \fdx{el_interpol_dow()@{\code{el\_interpol\_dow()}}} 753 %% 754 \kitem{dirichlet\_map(bound, bndry\_bits, mask)} 755 \fdx{dirichlet_map()@{\code{dirichlet\_map()}}}\hfill 756 757\end{descr} 758 759\subsubsection{BLAS-like Element-matrix and -vector operations} 760\label{S:elementblas} 761 762The source code listing below lists the proto-types, refer to 763\tableref{tab:elementblas1} and \tableref{tab:elementblas2} for a 764description of the respective operations. The routines in 765\tableref{tab:elementblas2} take an argument 766%% 767\bv 768\begin{lstlisting} 769const EL_SCHAR_VEC *bound. 770\end{lstlisting}\ev 771%% 772In this case the operations will act only on the rows $r$ which are 773not masked-out by \code{bound->vec[r] >= DIRICHLET}. The \code{bound} 774argument maybe \nil. 775%% 776\fdx{el_bi_mat_vec()@{\code{el\_bi\_mat\_vec()}}} 777\idx{BLAS for element vectors and matrices!el_bi_mat_vec()@{\code{el\_bi\_mat\_vec()}}} 778\fdx{el_bi_mat_vec_d()@{\code{el\_bi\_mat\_vec\_d()}}} 779\idx{BLAS for element vectors and matrices!el_bi_mat_vec_d()@{\code{el\_bi\_mat\_vec\_d()}}} 780\fdx{el_bi_mat_vec_dow()@{\code{el\_bi\_mat\_vec\_dow()}}} 781\idx{BLAS for element vectors and matrices!el_bi_mat_vec_dow()@{\code{el\_bi\_mat\_vec\_dow()}}} 782\fdx{el_bi_mat_vec_rrd()@{\code{el\_bi\_mat\_vec\_rrd()}}} 783\idx{BLAS for element vectors and matrices!el_bi_mat_vec_rrd()@{\code{el\_bi\_mat\_vec\_rrd()}}} 784\fdx{el_bi_mat_vec_scl_dow()@{\code{el\_bi\_mat\_vec\_scl\_dow()}}} 785\idx{BLAS for element vectors and matrices!el_bi_mat_vec_scl_dow()@{\code{el\_bi\_mat\_vec\_scl\_dow()}}} 786\fdx{el_bi_mat_vec_rdr()@{\code{el\_bi\_mat\_vec\_rdr()}}} 787\idx{BLAS for element vectors and matrices!el_bi_mat_vec_rdr()@{\code{el\_bi\_mat\_vec\_rdr()}}} 788\fdx{el_bi_mat_vec_dow_scl()@{\code{el\_bi\_mat\_vec\_dow\_scl()}}} 789\idx{BLAS for element vectors and matrices!el_bi_mat_vec_dow_scl()@{\code{el\_bi\_mat\_vec\_dow\_scl()}}} 790%% 791\fdx{el_gen_mat_vec()@{\code{el\_gen\_mat\_vec()}}} 792\idx{BLAS for element vectors and matrices!el_gen_mat_vec()@{\code{el\_gen\_mat\_vec()}}} 793\fdx{el_gen_mat_vec_d()@{\code{el\_gen\_mat\_vec\_d()}}} 794\idx{BLAS for element vectors and matrices!el_gen_mat_vec_d()@{\code{el\_gen\_mat\_vec\_d()}}} 795\fdx{el_gen_mat_vec_dow()@{\code{el\_gen\_mat\_vec\_dow()}}} 796\idx{BLAS for element vectors and matrices!el_gen_mat_vec_dow()@{\code{el\_gen\_mat\_vec\_dow()}}} 797\fdx{el_gen_mat_vec_rrd()@{\code{el\_gen\_mat\_vec\_rrd()}}} 798\idx{BLAS for element vectors and matrices!el_gen_mat_vec_rrd()@{\code{el\_gen\_mat\_vec\_rrd()}}} 799\fdx{el_gen_mat_vec_scl_dow()@{\code{el\_gen\_mat\_vec\_scl\_dow()}}} 800\idx{BLAS for element vectors and matrices!el_gen_mat_vec_scl_dow()@{\code{el\_gen\_mat\_vec\_scl\_dow()}}} 801\fdx{el_gen_mat_vec_rdr()@{\code{el\_gen\_mat\_vec\_rdr()}}} 802\idx{BLAS for element vectors and matrices!el_gen_mat_vec_rdr()@{\code{el\_gen\_mat\_vec\_rdr()}}} 803\fdx{el_gen_mat_vec_dow_scl()@{\code{el\_gen\_mat\_vec\_dow\_scl()}}} 804\idx{BLAS for element vectors and matrices!el_gen_mat_vec_dow_scl()@{\code{el\_gen\_mat\_vec\_dow\_scl()}}} 805%% 806\fdx{el_mat_vec()@{\code{el\_mat\_vec()}}} 807\idx{BLAS for element vectors and matrices!el_mat_vec()@{\code{el\_mat\_vec()}}} 808\fdx{el_mat_vec_d()@{\code{el\_mat\_vec\_d()}}} 809\idx{BLAS for element vectors and matrices!el_mat_vec_d()@{\code{el\_mat\_vec\_d()}}} 810\fdx{el_mat_vec_dow()@{\code{el\_mat\_vec\_dow()}}} 811\idx{BLAS for element vectors and matrices!el_mat_vec_dow()@{\code{el\_mat\_vec\_dow()}}} 812\fdx{el_mat_vec_rrd()@{\code{el\_mat\_vec\_rrd()}}} 813\idx{BLAS for element vectors and matrices!el_mat_vec_rrd()@{\code{el\_mat\_vec\_rrd()}}} 814\fdx{el_mat_vec_scl_dow()@{\code{el\_mat\_vec\_scl\_dow()}}} 815\idx{BLAS for element vectors and matrices!el_mat_vec_scl_dow()@{\code{el\_mat\_vec\_scl\_dow()}}} 816\fdx{el_mat_vec_rdr()@{\code{el\_mat\_vec\_rdr()}}} 817\idx{BLAS for element vectors and matrices!el_mat_vec_rdr()@{\code{el\_mat\_vec\_rdr()}}} 818\fdx{el_mat_vec_dow_scl()@{\code{el\_mat\_vec\_dow\_scl()}}} 819\idx{BLAS for element vectors and matrices!el_mat_vec_dow_scl()@{\code{el\_mat\_vec\_dow\_scl()}}} 820%% 821\fdx{bi_mat_el_vec()@{\code{bi\_mat\_el\_vec()}}} 822\idx{BLAS for element vectors and matrices!bi_mat_el_vec()@{\code{bi\_mat\_el\_vec()}}} 823\fdx{bi_mat_el_vec_d()@{\code{bi\_mat\_el\_vec\_d()}}} 824\idx{BLAS for element vectors and matrices!bi_mat_el_vec_d()@{\code{bi\_mat\_el\_vec\_d()}}} 825\fdx{bi_mat_el_vec_dow()@{\code{bi\_mat\_el\_vec\_dow()}}} 826\idx{BLAS for element vectors and matrices!bi_mat_el_vec_dow()@{\code{bi\_mat\_el\_vec\_dow()}}} 827\fdx{bi_mat_el_vec_rrd()@{\code{bi\_mat\_el\_vec\_rrd()}}} 828\idx{BLAS for element vectors and matrices!bi_mat_el_vec_rrd()@{\code{bi\_mat\_el\_vec\_rrd()}}} 829\fdx{bi_mat_el_vec_scl_dow()@{\code{bi\_mat\_el\_vec\_scl\_dow()}}} 830\idx{BLAS for element vectors and matrices!bi_mat_el_vec_scl_dow()@{\code{bi\_mat\_el\_vec\_scl\_dow()}}} 831\fdx{bi_mat_el_vec_rdr()@{\code{bi\_mat\_el\_vec\_rdr()}}} 832\idx{BLAS for element vectors and matrices!bi_mat_el_vec_rdr()@{\code{bi\_mat\_el\_vec\_rdr()}}} 833\fdx{bi_mat_el_vec_dow_scl()@{\code{bi\_mat\_el\_vec\_dow\_scl()}}} 834\idx{BLAS for element vectors and matrices!bi_mat_el_vec_dow_scl()@{\code{bi\_mat\_el\_vec\_dow\_scl()}}} 835%% 836\fdx{gen_mat_el_vec()@{\code{gen\_mat\_el\_vec()}}} 837\idx{BLAS for element vectors and matrices!gen_mat_el_vec()@{\code{gen\_mat\_el\_vec()}}} 838\fdx{gen_mat_el_vec_d()@{\code{gen\_mat\_el\_vec\_d()}}} 839\idx{BLAS for element vectors and matrices!gen_mat_el_vec_d()@{\code{gen\_mat\_el\_vec\_d()}}} 840\fdx{gen_mat_el_vec_dow()@{\code{gen\_mat\_el\_vec\_dow()}}} 841\idx{BLAS for element vectors and matrices!gen_mat_el_vec_dow()@{\code{gen\_mat\_el\_vec\_dow()}}} 842\fdx{gen_mat_el_vec_rrd()@{\code{gen\_mat\_el\_vec\_rrd()}}} 843\idx{BLAS for element vectors and matrices!gen_mat_el_vec_rrd()@{\code{gen\_mat\_el\_vec\_rrd()}}} 844\fdx{gen_mat_el_vec_scl_dow()@{\code{gen\_mat\_el\_vec\_scl\_dow()}}} 845\idx{BLAS for element vectors and matrices!gen_mat_el_vec_scl_dow()@{\code{gen\_mat\_el\_vec\_scl\_dow()}}} 846\fdx{gen_mat_el_vec_rdr()@{\code{gen\_mat\_el\_vec\_rdr()}}} 847\idx{BLAS for element vectors and matrices!gen_mat_el_vec_rdr()@{\code{gen\_mat\_el\_vec\_rdr()}}} 848\fdx{gen_mat_el_vec_dow_scl()@{\code{gen\_mat\_el\_vec\_dow\_scl()}}} 849\idx{BLAS for element vectors and matrices!gen_mat_el_vec_dow_scl()@{\code{gen\_mat\_el\_vec\_dow\_scl()}}} 850%% 851\fdx{mat_el_vec()@{\code{mat\_el\_vec()}}} 852\idx{BLAS for element vectors and matrices!mat_el_vec()@{\code{mat\_el\_vec()}}} 853\fdx{mat_el_vec_d()@{\code{mat\_el\_vec\_d()}}} 854\idx{BLAS for element vectors and matrices!mat_el_vec_d()@{\code{mat\_el\_vec\_d()}}} 855\fdx{mat_el_vec_dow()@{\code{mat\_el\_vec\_dow()}}} 856\idx{BLAS for element vectors and matrices!mat_el_vec_dow()@{\code{mat\_el\_vec\_dow()}}} 857\fdx{mat_el_vec_rrd()@{\code{mat\_el\_vec\_rrd()}}} 858\idx{BLAS for element vectors and matrices!mat_el_vec_rrd()@{\code{mat\_el\_vec\_rrd()}}} 859\fdx{mat_el_vec_scl_dow()@{\code{mat\_el\_vec\_scl\_dow()}}} 860\idx{BLAS for element vectors and matrices!mat_el_vec_scl_dow()@{\code{mat\_el\_vec\_scl\_dow()}}} 861\fdx{mat_el_vec_rdr()@{\code{mat\_el\_vec\_rdr()}}} 862\idx{BLAS for element vectors and matrices!mat_el_vec_rdr()@{\code{mat\_el\_vec\_rdr()}}} 863\fdx{mat_el_vec_dow_scl()@{\code{mat\_el\_vec\_dow\_scl()}}} 864\idx{BLAS for element vectors and matrices!mat_el_vec_dow_scl()@{\code{mat\_el\_vec\_dow\_scl()}}} 865%% 866\fdx{el_mat_set()@{\code{el\_mat\_set()}}} 867\idx{BLAS for element vectors and matrices!el_mat_set()@{\code{el\_mat\_set()}}} 868\fdx{el_mat_axey()@{\code{el\_mat\_axey()}}} 869\idx{BLAS for element vectors and matrices!el_mat_axey()@{\code{el\_mat\_axey()}}} 870\fdx{el_mat_axpy()@{\code{el\_mat\_axpy()}}} 871\idx{BLAS for element vectors and matrices!el_mat_axpy()@{\code{el\_mat\_axpy()}}} 872\fdx{el_mat_axpby()@{\code{el\_mat\_axpy()}}} 873\idx{BLAS for element vectors and matrices!el_mat_axpby()@{\code{el\_mat\_axpby()}}} 874\bv\begin{lstlisting} 875EL_REAL_VEC *el_bi_mat_vec(REAL a, const EL_MATRIX *A, 876 REAL b, const EL_MATRIX *B, 877 const EL_REAL_VEC *u_h, 878 REAL c, EL_REAL_VEC *f_h); 879EL_REAL_D_VEC *el_bi_mat_vec_d(REAL a, const EL_MATRIX *A, 880 REAL b, const EL_MATRIX *B, 881 const EL_REAL_D_VEC *u_h, 882 REAL c, EL_REAL_D_VEC *f_h); 883EL_REAL_VEC_D *el_bi_mat_vec_dow(REAL a, const EL_MATRIX *A, 884 REAL b, const EL_MATRIX *B, 885 const EL_REAL_VEC_D *u_h, 886 REAL c, EL_REAL_VEC_D *f_h); 887EL_REAL_VEC *el_bi_mat_vec_rrd(REAL a, const EL_MATRIX *A, 888 REAL b, const EL_MATRIX *B, 889 const EL_REAL_D_VEC *u_h, 890 REAL c, EL_REAL_VEC *f_h); 891EL_REAL_VEC *el_bi_mat_vec_scl_dow(REAL a, const EL_MATRIX *A, 892 REAL b, const EL_MATRIX *B, 893 const EL_REAL_VEC_D *u_h, 894 REAL c, EL_REAL_VEC *f_h); 895EL_REAL_D_VEC *el_bi_mat_vec_rdr(REAL a, const EL_MATRIX *A, 896 REAL b, const EL_MATRIX *B, 897 const EL_REAL_VEC *u_h, 898 REAL c, EL_REAL_D_VEC *f_h); 899EL_REAL_VEC_D *el_bi_mat_vec_dow_scl(REAL a, const EL_MATRIX *A, 900 REAL b, const EL_MATRIX *B, 901 const EL_REAL_VEC *u_h, 902 REAL c, EL_REAL_VEC_D *f_h); 903 904EL_REAL_VEC *el_gen_mat_vec(REAL a, const EL_MATRIX *A, 905 const EL_REAL_VEC *u_h, 906 REAL b, EL_REAL_VEC *f_h); 907EL_REAL_D_VEC *el_gen_mat_vec_d(REAL a, const EL_MATRIX *A, 908 const EL_REAL_D_VEC *u_h, 909 REAL b, EL_REAL_D_VEC *f_h); 910EL_REAL_VEC_D *el_gen_mat_vec_dow(REAL a, const EL_MATRIX *A, 911 const EL_REAL_VEC_D *u_h, 912 REAL b, EL_REAL_VEC_D *f_h); 913EL_REAL_VEC *el_gen_mat_vec_rrd(REAL a, const EL_MATRIX *A, 914 const EL_REAL_D_VEC *u_h, 915 REAL b, EL_REAL_VEC *f_h); 916EL_REAL_VEC *el_gen_mat_vec_scl_dow(REAL a, const EL_MATRIX *A, 917 const EL_REAL_VEC_D *u_h, 918 REAL b, EL_REAL_VEC *f_h); 919EL_REAL_D_VEC *el_gen_mat_vec_rdr(REAL a, const EL_MATRIX *A, 920 const EL_REAL_VEC *u_h, 921 REAL b, EL_REAL_D_VEC *f_h); 922EL_REAL_VEC_D *el_gen_mat_vec_dow_scl(REAL a, const EL_MATRIX *A, 923 const EL_REAL_VEC *u_h, 924 REAL b, EL_REAL_VEC_D *f_h); 925 926EL_REAL_VEC *el_mat_vec(REAL a, const EL_MATRIX *A, 927 const EL_REAL_VEC *u_h, EL_REAL_VEC *f_h); 928EL_REAL_D_VEC *el_mat_vec_d(REAL a, const EL_MATRIX *A, 929 const EL_REAL_D_VEC *u_h, EL_REAL_D_VEC *f_h); 930EL_REAL_VEC_D *el_mat_vec_dow(REAL a, const EL_MATRIX *A, 931 const EL_REAL_VEC_D *u_h, EL_REAL_VEC_D *f_h); 932EL_REAL_VEC *el_mat_vec_rrd(REAL a, const EL_MATRIX *A, 933 const EL_REAL_D_VEC *u_h, EL_REAL_VEC *f_h); 934EL_REAL_VEC *el_mat_vec_scl_dow(REAL a, const EL_MATRIX *A, 935 const EL_REAL_VEC_D *u_h, EL_REAL_VEC *f_h); 936EL_REAL_D_VEC *el_mat_vec_rdr(REAL a, const EL_MATRIX *A, 937 const EL_REAL_VEC *u_h, EL_REAL_D_VEC *f_h); 938EL_REAL_VEC_D *el_mat_vec_dow_scl(REAL a, const EL_MATRIX *A, 939 const EL_REAL_VEC *u_h, EL_REAL_VEC_D *f_h); 940 941void bi_mat_el_vec(REAL a, const EL_MATRIX *A, 942 REAL b, const EL_MATRIX *B, 943 const EL_REAL_VEC *u_h, REAL c, DOF_REAL_VEC *f_h, 944 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 945void bi_mat_el_vec_d(REAL a, const EL_MATRIX *A, 946 REAL b, const EL_MATRIX *B, 947 const EL_REAL_D_VEC *u_h, REAL c, DOF_REAL_D_VEC *f_h, 948 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 949void bi_mat_el_vec_dow(REAL a, const EL_MATRIX *A, 950 REAL b, const EL_MATRIX *B, 951 const EL_REAL_VEC_D *u_h, REAL c, DOF_REAL_VEC_D *f_h, 952 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 953void bi_mat_el_vec_rrd(REAL a, const EL_MATRIX *A, 954 REAL b, const EL_MATRIX *B, 955 const EL_REAL_D_VEC *u_h, REAL c, DOF_REAL_VEC *f_h, 956 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 957void bi_mat_el_vec_scl_dow(REAL a, const EL_MATRIX *A, 958 REAL b, const EL_MATRIX *B, 959 const EL_REAL_VEC_D *u_h, REAL c, DOF_REAL_VEC *f_h, 960 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 961void bi_mat_el_vec_rdr(REAL a, const EL_MATRIX *A, 962 REAL b, const EL_MATRIX *B, 963 const EL_REAL_VEC *u_h, REAL c, DOF_REAL_D_VEC *f_h, 964 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 965void bi_mat_el_vec_dow_scl(REAL a, const EL_MATRIX *A, 966 REAL b, const EL_MATRIX *B, 967 const EL_REAL_VEC *u_h, REAL c, DOF_REAL_VEC_D *f_h, 968 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 969 970void gen_mat_el_vec(REAL a, const EL_MATRIX *A, 971 const EL_REAL_VEC *u_h, REAL b, DOF_REAL_VEC *f_h, 972 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 973void gen_mat_el_vec_d(REAL a, const EL_MATRIX *A, 974 const EL_REAL_D_VEC *u_h, REAL b, DOF_REAL_D_VEC *f_h, 975 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 976void gen_mat_el_vec_dow(REAL a, const EL_MATRIX *A, 977 const EL_REAL_VEC_D *u_h, REAL b, DOF_REAL_VEC_D *f_h, 978 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 979void gen_mat_el_vec_rrd(REAL a, const EL_MATRIX *A, 980 const EL_REAL_D_VEC *u_h, REAL b, DOF_REAL_VEC *f_h, 981 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 982void gen_mat_el_vec_scl_dow(REAL a, const EL_MATRIX *A, 983 const EL_REAL_VEC_D *u_h, REAL b, DOF_REAL_VEC *f_h, 984 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 985void gen_mat_el_vec_rdr(REAL a, const EL_MATRIX *A, 986 const EL_REAL_VEC *u_h, REAL b, DOF_REAL_D_VEC *f_h, 987 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 988void gen_mat_el_vec_dow_scl(REAL a, const EL_MATRIX *A, 989 const EL_REAL_VEC *u_h, REAL b, DOF_REAL_VEC_D *f_h, 990 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 991 992void mat_el_vec(REAL a, const EL_MATRIX *A, 993 const EL_REAL_VEC *u_h, DOF_REAL_VEC *f_h, 994 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 995void mat_el_vec_d(REAL a, const EL_MATRIX *A, 996 const EL_REAL_D_VEC *u_h, DOF_REAL_D_VEC *f_h, 997 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 998void mat_el_vec_dow(REAL a, const EL_MATRIX *A, 999 const EL_REAL_VEC_D *u_h, DOF_REAL_VEC_D *f_h, 1000 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 1001void mat_el_vec_rrd(REAL a, const EL_MATRIX *A, 1002 const EL_REAL_D_VEC *u_h, DOF_REAL_VEC *f_h, 1003 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 1004void mat_el_vec_scl_dow(REAL a, const EL_MATRIX *A, 1005 const EL_REAL_VEC_D *u_h, DOF_REAL_VEC *f_h, 1006 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 1007void mat_el_vec_rdr(REAL a, const EL_MATRIX *A, 1008 const EL_REAL_VEC *u_h, DOF_REAL_D_VEC *f_h, 1009 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 1010void mat_el_vec_dow_scl(REAL a, const EL_MATRIX *A, 1011 const EL_REAL_VEC *u_h, DOF_REAL_VEC_D *f_h, 1012 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound); 1013 1014EL_MATRIX *el_mat_set(REAL a, EL_MATRIX *result); 1015EL_MATRIX *el_mat_axey(REAL a, const EL_MATRIX *A, EL_MATRIX *result); 1016EL_MATRIX *el_mat_axpy(REAL a, const EL_MATRIX *A, EL_MATRIX *result); 1017EL_MATRIX *el_mat_axpby(REAL a, const EL_MATRIX *A, 1018 REAL b, const EL_MATRIX *B, EL_MATRIX *result); 1019\end{lstlisting}\ev 1020 1021\begin{table}[htbp] 1022\idx{BLAS for element vectors and matrices} 1023\fdx{el_bi_mat_vec()@{\code{el\_bi\_mat\_vec()}}} 1024\idx{BLAS for element vectors and matrices!el_bi_mat_vec()@{\code{el\_bi\_mat\_vec()}}} 1025\fdx{el_bi_mat_vec_d()@{\code{el\_bi\_mat\_vec\_d()}}} 1026\idx{BLAS for element vectors and matrices!el_bi_mat_vec_d()@{\code{el\_bi\_mat\_vec\_d()}}} 1027\fdx{el_bi_mat_vec_dow()@{\code{el\_bi\_mat\_vec\_dow()}}} 1028\idx{BLAS for element vectors and matrices!el_bi_mat_vec_dow()@{\code{el\_bi\_mat\_vec\_dow()}}} 1029\fdx{el_bi_mat_vec_rrd()@{\code{el\_bi\_mat\_vec\_rrd()}}} 1030\idx{BLAS for element vectors and matrices!el_bi_mat_vec_rrd()@{\code{el\_bi\_mat\_vec\_rrd()}}} 1031\fdx{el_bi_mat_vec_scl_dow()@{\code{el\_bi\_mat\_vec\_scl\_dow()}}} 1032\idx{BLAS for element vectors and matrices!el_bi_mat_vec_scl_dow()@{\code{el\_bi\_mat\_vec\_scl\_dow()}}} 1033\fdx{el_bi_mat_vec_rdr()@{\code{el\_bi\_mat\_vec\_rdr()}}} 1034\idx{BLAS for element vectors and matrices!el_bi_mat_vec_rdr()@{\code{el\_bi\_mat\_vec\_rdr()}}} 1035\fdx{el_bi_mat_vec_dow_scl()@{\code{el\_bi\_mat\_vec\_dow\_scl()}}} 1036\idx{BLAS for element vectors and matrices!el_bi_mat_vec_dow_scl()@{\code{el\_bi\_mat\_vec\_dow\_scl()}}} 1037%% 1038\fdx{el_gen_mat_vec()@{\code{el\_gen\_mat\_vec()}}} 1039\idx{BLAS for element vectors and matrices!el_gen_mat_vec()@{\code{el\_gen\_mat\_vec()}}} 1040\fdx{el_gen_mat_vec_d()@{\code{el\_gen\_mat\_vec\_d()}}} 1041\idx{BLAS for element vectors and matrices!el_gen_mat_vec_d()@{\code{el\_gen\_mat\_vec\_d()}}} 1042\fdx{el_gen_mat_vec_dow()@{\code{el\_gen\_mat\_vec\_dow()}}} 1043\idx{BLAS for element vectors and matrices!el_gen_mat_vec_dow()@{\code{el\_gen\_mat\_vec\_dow()}}} 1044\fdx{el_gen_mat_vec_rrd()@{\code{el\_gen\_mat\_vec\_rrd()}}} 1045\idx{BLAS for element vectors and matrices!el_gen_mat_vec_rrd()@{\code{el\_gen\_mat\_vec\_rrd()}}} 1046\fdx{el_gen_mat_vec_scl_dow()@{\code{el\_gen\_mat\_vec\_scl\_dow()}}} 1047\idx{BLAS for element vectors and matrices!el_gen_mat_vec_scl_dow()@{\code{el\_gen\_mat\_vec\_scl\_dow()}}} 1048\fdx{el_gen_mat_vec_rdr()@{\code{el\_gen\_mat\_vec\_rdr()}}} 1049\idx{BLAS for element vectors and matrices!el_gen_mat_vec_rdr()@{\code{el\_gen\_mat\_vec\_rdr()}}} 1050\fdx{el_gen_mat_vec_dow_scl()@{\code{el\_gen\_mat\_vec\_dow\_scl()}}} 1051\idx{BLAS for element vectors and matrices!el_gen_mat_vec_dow_scl()@{\code{el\_gen\_mat\_vec\_dow\_scl()}}} 1052%% 1053\fdx{el_mat_vec()@{\code{el\_mat\_vec()}}} 1054\idx{BLAS for element vectors and matrices!el_mat_vec()@{\code{el\_mat\_vec()}}} 1055\fdx{el_mat_vec_d()@{\code{el\_mat\_vec\_d()}}} 1056\idx{BLAS for element vectors and matrices!el_mat_vec_d()@{\code{el\_mat\_vec\_d()}}} 1057\fdx{el_mat_vec_dow()@{\code{el\_mat\_vec\_dow()}}} 1058\idx{BLAS for element vectors and matrices!el_mat_vec_dow()@{\code{el\_mat\_vec\_dow()}}} 1059\fdx{el_mat_vec_rrd()@{\code{el\_mat\_vec\_rrd()}}} 1060\idx{BLAS for element vectors and matrices!el_mat_vec_rrd()@{\code{el\_mat\_vec\_rrd()}}} 1061\fdx{el_mat_vec_scl_dow()@{\code{el\_mat\_vec\_scl\_dow()}}} 1062\idx{BLAS for element vectors and matrices!el_mat_vec_scl_dow()@{\code{el\_mat\_vec\_scl\_dow()}}} 1063\fdx{el_mat_vec_rdr()@{\code{el\_mat\_vec\_rdr()}}} 1064\idx{BLAS for element vectors and matrices!el_mat_vec_rdr()@{\code{el\_mat\_vec\_rdr()}}} 1065\fdx{el_mat_vec_dow_scl()@{\code{el\_mat\_vec\_dow\_scl()}}} 1066\idx{BLAS for element vectors and matrices!el_mat_vec_dow_scl()@{\code{el\_mat\_vec\_dow\_scl()}}} 1067%% 1068\fdx{el_mat_set()@{\code{el\_mat\_set()}}} 1069\idx{BLAS for element vectors and matrices!el_mat_set()@{\code{el\_mat\_set()}}} 1070\fdx{el_mat_axey()@{\code{el\_mat\_axey()}}} 1071\idx{BLAS for element vectors and matrices!el_mat_axey()@{\code{el\_mat\_axey()}}} 1072\fdx{el_mat_axpy()@{\code{el\_mat\_axpy()}}} 1073\idx{BLAS for element vectors and matrices!el_mat_axpy()@{\code{el\_mat\_axpy()}}} 1074\fdx{el_mat_axpby()@{\code{el\_mat\_axpby()}}} 1075\idx{BLAS for element vectors and matrices!el_mat_axpby()@{\code{el\_mat\_axpby()}}} 1076\begin{center}{\small 1077 \begin{tabular}{|l|l|} 1078 \hline 1079 \Strut\verb|f = el_mat_vec(a, A, u, f)| & 1080 $f_i \leftarrow (a\,A\,u)_i$ \\ 1081 \Strut\verb|f = el_mat_vec_d(a, A, u, f)| & \\ 1082 \Strut\verb|f = el_mat_vec_dow(a, A, u, f)| & \\ 1083 \Strut\verb|f = el_mat_vec_rrd(a, A, u, f)| & \\ 1084 \Strut\verb|f = el_mat_vec_scl_dow(a, A, u, f)| &\\ 1085 \Strut\verb|f = el_mat_vec_rdr(a, A, u, f)| &\\ 1086 \Strut\verb|f = el_mat_vec_dow_scl(a, A, u, f)|&\\ 1087 \hline 1088 \Strut\verb|f = el_gen_mat_vec(a, A, u, b, f)| & 1089 $f_i \leftarrow (a\,A\,u + b\,f)_i$ \\ 1090 \Strut\verb|f = el_gen_mat_vec_d(a, A, u, b, f)| & \\ 1091 \Strut\verb|f = el_gen_mat_vec_dow(a, A, u, b, f)| & \\ 1092 \Strut\verb|f = el_gen_mat_vec_rrd(a, A, u, b, f)| & \\ 1093 \Strut\verb|f = el_gen_mat_vec_scl_dow(a, A, u, b, f)| &\\ 1094 \Strut\verb|f = el_gen_mat_vec_rdr(a, A, u, b, f)| &\\ 1095 \Strut\verb|f = el_gen_mat_vec_dow_scl(a, A, u, b, f)|&\\ 1096 \hline 1097 \Strut\verb|f = el_bi_mat_vec(a, A, b, B, u, c, f)| & 1098 $f_i \leftarrow ((a\,A + b\,B)\,u + c\,f)_i$ \\ 1099 \Strut\verb|f = el_bi_mat_vec_d(a, A, b, B, u, c, f)| & \\ 1100 \Strut\verb|f = el_bi_mat_vec_dow(a, A, b, B, u, c, f)| & \\ 1101 \Strut\verb|f = el_bi_mat_vec_rrd(a, A, b, B, u, c, f)| & \\ 1102 \Strut\verb|f = el_bi_mat_vec_scl_dow(a, A, b, B, u, c, f)| &\\ 1103 \Strut\verb|f = el_bi_mat_vec_rdr(a, A, b, B, u, c, f)| &\\ 1104 \Strut\verb|f = el_bi_mat_vec_dow_scl(a, A, b, B, u, c, f)|&\\ 1105 \hline 1106 \Strut\verb|A = el_mat_set(a, A)| & 1107 $A_{ij} \leftarrow a$\\ 1108 \Strut\verb|B = el_mat_axey(a, A, B)| & 1109 $B_{ij} \leftarrow a\,A_{ij}$\\ 1110 \Strut\verb|B = el_mat_axpy(a, A, B)| & 1111 $B_{ij} \leftarrow a\,A_{ij}+B_{ij}$\\ 1112 \Strut\verb|C = el_mat_axpby(a, A, b, B, C)| & 1113 $C_{ij} \leftarrow a\,A_{ij}+b\,B_{ij}$\\ 1114 \hline 1115 \end{tabular}}\end{center} 1116\caption[BLAS-operations for element-vectors and -matrices] 1117{BLAS-operations for element-vectors and -matrices. $A$ and $B$ denote element matrices, $u$ and $f$ element vectors, $a$, $b$, $c$ are numbers.} 1118\label{tab:elementblas1} 1119\end{table} 1120 1121\begin{table}[htbp] 1122\idx{BLAS for element vectors and matrices} 1123\fdx{bi_mat_el_vec()@{\code{bi\_mat\_el\_vec()}}} 1124\idx{BLAS for element vectors and matrices!bi_mat_el_vec()@{\code{bi\_mat\_el\_vec()}}} 1125\fdx{bi_mat_el_vec_d()@{\code{bi\_mat\_el\_vec\_d()}}} 1126\idx{BLAS for element vectors and matrices!bi_mat_el_vec_d()@{\code{bi\_mat\_el\_vec\_d()}}} 1127\fdx{bi_mat_el_vec_dow()@{\code{bi\_mat\_el\_vec\_dow()}}} 1128\idx{BLAS for element vectors and matrices!bi_mat_el_vec_dow()@{\code{bi\_mat\_el\_vec\_dow()}}} 1129\fdx{bi_mat_el_vec_rrd()@{\code{bi\_mat\_el\_vec\_rrd()}}} 1130\idx{BLAS for element vectors and matrices!bi_mat_el_vec_rrd()@{\code{bi\_mat\_el\_vec\_rrd()}}} 1131\fdx{bi_mat_el_vec_scl_dow()@{\code{bi\_mat\_el\_vec\_scl\_dow()}}} 1132\idx{BLAS for element vectors and matrices!bi_mat_el_vec_scl_dow()@{\code{bi\_mat\_el\_vec\_scl\_dow()}}} 1133\fdx{bi_mat_el_vec_rdr()@{\code{bi\_mat\_el\_vec\_rdr()}}} 1134\idx{BLAS for element vectors and matrices!bi_mat_el_vec_rdr()@{\code{bi\_mat\_el\_vec\_rdr()}}} 1135\fdx{bi_mat_el_vec_dow_scl()@{\code{bi\_mat\_el\_vec\_dow\_scl()}}} 1136\idx{BLAS for element vectors and matrices!bi_mat_el_vec_dow_scl()@{\code{bi\_mat\_el\_vec\_dow\_scl()}}} 1137%% 1138\fdx{gen_mat_el_vec()@{\code{gen\_mat\_el\_vec()}}} 1139\idx{BLAS for element vectors and matrices!gen_mat_el_vec()@{\code{gen\_mat\_el\_vec()}}} 1140\fdx{gen_mat_el_vec_d()@{\code{gen\_mat\_el\_vec\_d()}}} 1141\idx{BLAS for element vectors and matrices!gen_mat_el_vec_d()@{\code{gen\_mat\_el\_vec\_d()}}} 1142\fdx{gen_mat_el_vec_dow()@{\code{gen\_mat\_el\_vec\_dow()}}} 1143\idx{BLAS for element vectors and matrices!gen_mat_el_vec_dow()@{\code{gen\_mat\_el\_vec\_dow()}}} 1144\fdx{gen_mat_el_vec_rrd()@{\code{gen\_mat\_el\_vec\_rrd()}}} 1145\idx{BLAS for element vectors and matrices!gen_mat_el_vec_rrd()@{\code{gen\_mat\_el\_vec\_rrd()}}} 1146\fdx{gen_mat_el_vec_scl_dow()@{\code{gen\_mat\_el\_vec\_scl\_dow()}}} 1147\idx{BLAS for element vectors and matrices!gen_mat_el_vec_scl_dow()@{\code{gen\_mat\_el\_vec\_scl\_dow()}}} 1148\fdx{gen_mat_el_vec_rdr()@{\code{gen\_mat\_el\_vec\_rdr()}}} 1149\idx{BLAS for element vectors and matrices!gen_mat_el_vec_rdr()@{\code{gen\_mat\_el\_vec\_rdr()}}} 1150\fdx{gen_mat_el_vec_dow_scl()@{\code{gen\_mat\_el\_vec\_dow\_scl()}}} 1151\idx{BLAS for element vectors and matrices!gen_mat_el_vec_dow_scl()@{\code{gen\_mat\_el\_vec\_dow\_scl()}}} 1152%% 1153\fdx{mat_el_vec()@{\code{mat\_el\_vec()}}} 1154\idx{BLAS for element vectors and matrices!mat_el_vec()@{\code{mat\_el\_vec()}}} 1155\fdx{mat_el_vec_d()@{\code{mat\_el\_vec\_d()}}} 1156\idx{BLAS for element vectors and matrices!mat_el_vec_d()@{\code{mat\_el\_vec\_d()}}} 1157\fdx{mat_el_vec_dow()@{\code{mat\_el\_vec\_dow()}}} 1158\idx{BLAS for element vectors and matrices!mat_el_vec_dow()@{\code{mat\_el\_vec\_dow()}}} 1159\fdx{mat_el_vec_rrd()@{\code{mat\_el\_vec\_rrd()}}} 1160\idx{BLAS for element vectors and matrices!mat_el_vec_rrd()@{\code{mat\_el\_vec\_rrd()}}} 1161\fdx{mat_el_vec_scl_dow()@{\code{mat\_el\_vec\_scl\_dow()}}} 1162\idx{BLAS for element vectors and matrices!mat_el_vec_scl_dow()@{\code{mat\_el\_vec\_scl\_dow()}}} 1163\fdx{mat_el_vec_rdr()@{\code{mat\_el\_vec\_rdr()}}} 1164\idx{BLAS for element vectors and matrices!mat_el_vec_rdr()@{\code{mat\_el\_vec\_rdr()}}} 1165\fdx{mat_el_vec_dow_scl()@{\code{mat\_el\_vec\_dow\_scl()}}} 1166\idx{BLAS for element vectors and matrices!mat_el_vec_dow_scl()@{\code{mat\_el\_vec\_dow\_scl()}}} 1167\begin{center}{\small 1168 \begin{tabular}{|l|l|} 1169 \hline 1170 \Strut\verb|f = mat_el_vec(a, A, u, f, dof, mask)| & 1171 $f[dof[i]] \leftarrow (a\,A\,u)_i$ \\ 1172 \Strut\verb|mat_el_vec_d(a, A, u, f, dof, mask)| & 1173 if \code{mask[i] != DIRICHLET}\\ 1174 \Strut\verb|mat_el_vec_dow(a, A, u, f, dof, mask)| & 1175 or \code{mask == \nil}\\ 1176 \Strut\verb|mat_el_vec_rrd(a, A, u, f, dof, mask)| & \\ 1177 \Strut\verb|mat_el_vec_scl_dow(a, A, u, f, dof, mask)| &\\ 1178 \Strut\verb|mat_el_vec_rdr(a, A, u, f, dof, mask)| &\\ 1179 \Strut\verb|mat_el_vec_dow_scl(a, A, u, f, dof, mask)|&\\ 1180 \hline 1181 \Strut\verb|gen_mat_el_vec(a, A, u, b, f, dof, mask)| & 1182 $f[dof[i]] \leftarrow (a\,A\,u)_i + b\,f[dof[i]]$ \\ 1183 \Strut\verb|gen_mat_el_vec_d(a, A, u, b, f, dof, mask)| & 1184 if \code{mask[i] != DIRICHLET}\\ 1185 \Strut\verb|gen_mat_el_vec_dow(a, A, u, b, f, dof, mask)| & 1186 or \code{mask == \nil}\\ 1187 \Strut\verb|gen_mat_el_vec_rrd(a, A, u, b, f, dof, mask)| & \\ 1188 \Strut\verb|gen_mat_el_vec_scl_dow(a, A, u, b, f, dof, mask)| &\\ 1189 \Strut\verb|gen_mat_el_vec_rdr(a, A, u, b, f, dof, mask)| &\\ 1190 \Strut\verb|gen_mat_el_vec_dow_scl(a, A, u, b, f, dof, mask)|&\\ 1191 \hline 1192 \Strut\verb|bi_mat_el_vec(a, A, b, B, u, c, f, dof, mask)| & 1193 $f[dof[i]]$ \\ 1194 \Strut\verb|bi_mat_el_vec_d(a, A, b, B, u, c, f, dof, mask)| & 1195 $\leftarrow ((a\,A + b\,B)\,u)_i + c\,f[dof[i]]$ \\ 1196 \Strut\verb|bi_mat_el_vec_dow(a, A, b, B, u, c, f, dof, mask)| & 1197 if \code{mask[i] != DIRICHLET}\\ 1198 \Strut\verb|bi_mat_el_vec_rrd(a, A, b, B, u, c, f, dof, mask)| & 1199 or \code{mask == \nil}\\ 1200 \Strut\verb|bi_mat_el_vec_scl_dow(a, A, b, B, u, c, f, dof, mask)| &\\ 1201 \Strut\verb|bi_mat_el_vec_rdr(a, A, b, B, u, c, f, dof, mask)| &\\ 1202 \Strut\verb|bi_mat_el_vec_dow_scl(a, A, b, B, u, c, f, dof, mask)|&\\ 1203 \hline 1204\end{tabular}}\end{center} 1205\caption[BLAS-operations for element-vectors and -matrices] 1206{BLAS-operations for element-vectors and -matrices. $A$ and $B$ denote element matrices, $u$ an element vector. $f$ is a global \code{DOF}-vector. \code{mask} is an \code{EL\_SCHAR\_VEC} masking out certain \emph{local} \code{DOF}s. \code{mask} may be \nil. \code{dof} is \code{EL\_DOF\_VEC} mapping local to global \code{DOF}s. $a$, $b$, $c$ are numbers.} 1207\label{tab:elementblas2} 1208\end{table} 1209 1210\subsection{Data structures and functions for matrix assemblage}% 1211\label{S:matrix_assemblage}% 1212 1213The following structure holds full information for the assembling 1214of scalar element matrices. This structure is used by the function 1215\code{update\_matrix()} described below. 1216\ddx{EL_MATRIX_INFO@{\code{EL\_MATRIX\_INFO}}} 1217\idx{assemblage tools!EL_MATRIX_INFO@{\code{EL\_MATRIX\_INFO}}} 1218%% 1219\newcommand{\ELMATRIXINFO}{\hyperref[T:EL_MATRIX_INFO]{\code{EL\_MATRIX\_INFO}}\xspace} 1220%% 1221\bv\begin{lstlisting}[label=T:EL_MATRIX_INFO,name=EL_MATRIX_INFO] 1222typedef struct el_matrix_info EL_MATRIX_INFO; 1223 1224struct el_matrix_info 1225{ 1226 const FE_SPACE *row_fe_space; 1227 const FE_SPACE *col_fe_space; 1228 1229 MATENT_TYPE krn_blk_type; 1230 1231 BNDRY_FLAGS dirichlet_bndry; 1232 REAL factor; 1233 1234 EL_MATRIX_FCT el_matrix_fct; 1235 void *fill_info; 1236 1237 const EL_MATRIX_FCT *neigh_el_mat_fcts; 1238 void *neigh_fill_info; 1239 1240 FLAGS fill_flag; 1241}; 1242\end{lstlisting}\ev 1243Description: 1244\begin{descr} 1245 \kitem{row\_fe\_space} pointer to a finite element space connected 1246 to the row DOFs of the resulting matrix. 1247 %% 1248 \kitem{col\_fe\_space} pointer to a finite element space connected 1249 to the columns DOFs of the resulting matrix. 1250 %% 1251 \kitem{krn\_blk\_type} defines the block-matrix type of matrix entries 1252 %% 1253 \kitem{dirichlet\_boundary} bndry-type bit-mask for Dirichlet-boundary 1254 conditions built into the matrix 1255 %% 1256 \kitem{factor} is a multiplier for the element contributions; usually 1257 \code{factor} is \code{1} or \code{-1}. 1258 %% 1259 \kitem{el\_matrix\_fct}is a pointer to a function for the computation 1260 of the element matrix; \code{el\_matrix\_fct(el\_info, fill\_info)} 1261 returns a pointer to a matrix of size 1262 $\code{n\_row}\times\code{n\_col}$ storing the element matrix 1263 on element \code{el\_info->el}; \code{fill\_info} is a pointer 1264 to data needed by \code{el\_matrix\_fct()}; the function has 1265 to provide memory for storing the element matrix, which 1266 can be overwritten on the next call. 1267 %% 1268 \kitem{fill\_info} pointer to data needed by \code{el\_matrix\_fct()}; 1269 will be given as second argument to this function. 1270 %% 1271 \kitem{neigh\_el\_mat\_fcts} If the \code{BNDRY\_OPERATOR\_INFO} 1272 (code-listing \ref{T:BNDRY_OPERATOR_INFO}) structure passed to 1273 \code{fill\_matrix\_info()} was flagged as discontinuous, then this 1274 is the base-address of an array storing \code{N\_NEIGH(mesh->dim)} 1275 many element-matrix functions which pair the local basis-function 1276 set with the local basis function set on the neighbor. 1277 Intentionally, this is meant to support assembling linear systems in 1278 the context of DG-methods. The idea is that 1279 %% 1280 \code{EL\_MATRIX\_INFO.neigh\_el\_mat\_fcts[neigh\_nr](el\_info, EL\_MATRIX\_INFO.neigh\_fill\_info)} 1281 %% 1282 assembles a jump-term where the local basis functions on the element 1283 described by \code{el\_info} are used as test-functions 1284 (corresponding to the rows of the element matrix) and the local 1285 basis function set on the neighbour element defines the local space 1286 of ansatz-functions (column-space). 1287 %% 1288 \kitem{neigh\_fill\_info} Data pointer passed to the element-matrix 1289 functions stored in \code{neigh\_el\_mat\_fcts}. 1290 %% 1291 \kitem{fill\_flag}the flag for the mesh traversal for assembling the 1292 matrix. 1293\end{descr} 1294The following function updates a matrix by assembling element 1295contributions during mesh traversal; information for computing 1296the element matrices is provided in an \code{EL\_MATRIX\_INFO} structure: 1297\fdx{update_matrix()@{\code{update\_matrix()}}} 1298\idx{assemblage tools!update_matrix()@{\code{update\_matrix(}}} 1299\bv\begin{lstlisting}[name=update_matrix(),label=T:update_matrix_proto] 1300void update_matrix(DOF_MATRIX *dof_matrix, const EL_MATRIX_INFO *minfo, 1301 MatrixTranspose transpose);; 1302\end{lstlisting}\ev 1303Description: 1304\begin{descr} 1305 \kitem{update\_matrix(matrix, info, transpose)} updates the matrix 1306 \code{matrix} by traversing the underlying mesh and assembling the 1307 element contributions into the matrix; information about the 1308 computation of element matrices and connection of local and global 1309 DOFs is stored in \code{info}. 1310 1311 The flags for the mesh traversal of the mesh 1312 \code{matrix->fe\_space->mesh} are stored at \code{info->fill\_flag} 1313 which specifies the elements to be visited and information that 1314 should be present on the elements for the calculation of the element 1315 matrices and boundary information (if \code{info->get\_bound} is not 1316 \nil). 1317 1318 On the elements, information about the row DOFs is accessed by 1319 \code{info->get\_row\_dof} using \code{info->row\_admin}; this 1320 vector is also used for the column DOFs if \code{info->n\_col} is 1321 less or equal to zero, or \code{info->get\_col\_admin} or 1322 \code{info->get\_col\_dof} is a \nil pointer; when row and column 1323 DOFs are the same, the boundary type of the DOFs is accessed by 1324 \code{info->get\_bound} if \code{info->get\_bound} is not a \nil 1325 pointer; then the element matrix is computed by 1326 \code{info->el\_matrix\_fct(el\_info, info->fill\_info)}; these 1327 contributions, multiplied by \code{info->factor}, are eventually 1328 added to \code{matrix} by a call of \code{add\_element\_matrix()} 1329 with all information about row and column DOFs, the element matrix, 1330 and boundary types, if available. 1331 1332 \code{update\_matrix()} acts \emph{additive}, the 1333 element-contributions are added to the data already present in 1334 \code{dof\_matrix}. This makes several calls for the assemblage of 1335 one matrix possible. \code{clear\_dof\_matrix()} can be used to 1336 erase the contents of \code{dof\_matrix} prior to calling 1337 \code{update\_matrix()}. 1338\begin{description} 1339\item[Parameters]\hfill 1340 \begin{descr} 1341 \kitem{dof\_matrix} The global \code{DOF\_MATRIX} to add data to. 1342 %% 1343 \kitem{minfo} The element-matrix handle, as returned by 1344 \code{fill\_matrix\_info()} or 1345 %% 1346 \kitem{transpose} 1347 \end{descr} 1348\end{description} 1349\end{descr} 1350 1351\subsection{Matrix assemblage for second order problems}% 1352\label{S:matrix_assemblage_scalar}% 1353 1354Now we want to describe some tools which enable an easy assemblage of 1355the system matrix in the case of scalar elliptic problems. For this we 1356have to provide a function for the calculation of the element 1357matrix. For a general scalar problem the element matrix 1358$\vL_S = (L_S^{ij})_{i,j=1,\dots,m}$ is given by (recall 1359\mathref{book:E:L-phii-phij} on page 1360\pageref{book:E:L-phii-phij}) 1361%in \secref{book:S:Dis2ndorder} 1362\begin{align*} 1363L_S^{ij} &= 1364\int_{\Shat} 1365 \nablal \pbar^{i}(\lambda(\xhat)) \cdot \bar A(\lambda(\xhat))\, 1366 \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat 1367+ 1368\int_{\Shat} \pbar^{i}(\lambda(\xhat))\; 1369 \bar b(\lambda(\xhat)) \cdot 1370 \nablal \pbar^{j}(\lambda(\xhat)) \,d\xhat\\ 1371&\qquad + 1372\int_{\Shat} \bar c(\lambda(\xhat))\, \pbar^{i}(\lambda(\xhat)) \, 1373 \pbar^{j}(\lambda(\xhat))\,d\xhat, 1374\end{align*} 1375where $\bar A$, $\bar b$, and $\bar c$ are functions depending 1376on given data and on the actual element, namely 1377\begin{align*} 1378\bar A(\lambda) &:= 1379\left(\bar a_{kl}(\lambda)\right)_{k,l = 0,\ldots,d} 1380 := |\det DF_S(\xhat(\lambda))| \, 1381\Lambda(x(\lambda)) \, A(x(\lambda)) \, \Lambda^t(x(\lambda)),\\ 1382\bar b(\lambda) & := 1383\left(\bar b_{l}(\lambda)\right)_{l = 0,\ldots,d} 1384 := |\det DF_S(\xhat(\lambda))| \, 1385\Lambda(x(\lambda)) \, b(x(\lambda)), \quad \mbox{and}\\ 1386\bar c(\lambda) & := |\det DF_S(\xhat(\lambda))| \, c(x(\lambda)). 1387\end{align*} 1388Having access to functions for the evaluation of $\bar A$, $\bar b$, 1389and $\bar c$ at given quadrature nodes, the above integrals can be 1390computed by some general routine for any set of local basis functions 1391using quadrature. Additionally, if a coefficient is piecewise constant 1392on the mesh, only an integration of basis functions has to be done 1393(compare \mathref{book:E:quad_LS} on page \pageref{book:E:quad_LS}) 1394for this term. Here we can use pre--computed integrals of the basis 1395functions on the standard element and transform th-em to the actual 1396element. Such a computation is usually much faster than using 1397quadrature on each single element. Data structures for storing such 1398pre--computed values are described in \secref{S:ass_info}. 1399 1400For the assemblage routines which we will describe now, we use 1401the following slight generalization: In the discretization of 1402the first order term, sometimes integration by parts is used 1403too. For a divergence free vector field 1404$b$ and purely Dirichlet boundary values this leads for instance to 1405\[ 1406\int_\Omega \varphi(x)\, b(x) \cdot \nabla u(x)\,dx = 1407\frac12 \left(\int_\Omega \varphi(x)\, b(x) \cdot \nabla u(x)\,dx 1408- \int_\Omega \nabla \varphi(x) \cdot b(x) \, u(x)\,dx\right) 1409\] 1410yielding a modified first order term for the element matrix 1411\[ 1412\int_{\Shat} \pbar^{i}(\lambda(\xhat))\; 1413 \frac12 \bar b(\lambda(\xhat)) \cdot 1414 \nablal \pbar^{j}(\lambda(\xhat)) \,d\xhat 1415- \int_{\Shat} \nablal \pbar^{i}(\lambda(\xhat))\cdot 1416 \frac12 \bar b(\lambda(\xhat)) \; \pbar^{j}(\lambda(\xhat)) \,d\xhat. 1417\] 1418Secondly, we allow that we have two finite element spaces with local 1419basis functions $\{\bar\psi_i\}_{i=1,\dots,n}$ and 1420$\{\bar\vphi_i\}_{i=1,\dots,m}$. 1421 1422In general the following contributions of the element matrix 1423$\vL_S=(L_S^{ij})_{\substack{i=1,\dots,n\\j=1,\dots,m}}$ have to be 1424computed: 1425\begin{align*} 1426&\int_{\Shat} 1427 \nablal \bar\psi^{i}(\lambda(\xhat)) \cdot \bar A(\lambda(\xhat))\, 1428 \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat & &\mbox{second order term,}\\ 1429&\begin{array}{l}% 1430\ds \int_{\Shat} \bar\psi^{i}(\lambda(\xhat))\; 1431 \bar b^0(\lambda(\xhat)) \cdot \nablal \pbar^{j}(\lambda(\xhat)) \,d\xhat 1432\\[3mm] 1433\ds \int_{\Shat}^{} \nablal \bar\psi^{i}(\lambda(\xhat)) \cdot\; 1434 \bar b^1(\lambda(\xhat))\, \pbar^{j}(\lambda(\xhat)) \,d\xhat\\ 1435\end{array} 1436& &\mbox{first order terms,}\\ 1437& 1438\int_{\Shat} \bar c(\lambda(\xhat))\, \bar\psi^{i}(\lambda(\xhat)) \, 1439 \pbar^{j}(\lambda(\xhat))\,d\xhat 1440& &\mbox{zero order term,} 1441\end{align*} 1442where for instance $\bar b^0 = \bar b$ and $\bar b^1 = 0$, or using 1443integration by parts $\bar b^0 = \frac12\bar b$ and 1444$\bar b^1=-\frac12\bar b$. 1445 1446In order to store information about the finite element spaces, the 1447problem dependent functions $\bar A$, $\bar b^0$, $\bar b^1$, $\bar c$ 1448and the quadrature that should be used for the numerical 1449integration of the element matrix, we define the following data 1450structure: 1451%% 1452\ddx{OPERATOR_INFO@{\code{OPERATOR\_INFO}}} 1453\idx{assemblage tools!OPERATOR_INFO@{\code{OPERATOR\_INFO}}} 1454%% 1455\newcommand{\OPERATORINFO}{\hyperref[T:OPERATOR_INFO]{\code{OPERATOR\_INFO}}\xspace} 1456%% 1457\bv\begin{lstlisting}[label=T:OPERATOR_INFO,name=OPERATOR_INFO] 1458typedef struct operator_info OPERATOR_INFO; 1459 1460struct operator_info 1461{ 1462 const FE_SPACE *row_fe_space; /* range fe-space */ 1463 const FE_SPACE *col_fe_space; /* domain fe-space */ 1464 1465 const QUAD *quad[3]; 1466 1467 bool (*init_element)(const EL_INFO *el_info, 1468 const QUAD *quad[3], void *apd); 1469 1470 union { 1471 const REAL_B *(*real)(const EL_INFO *el_info, 1472 const QUAD *quad, int iq, void *apd); 1473 const REAL_BD *(*real_d)(const EL_INFO *el_info, 1474 const QUAD *quad, int iq, void *apd); 1475 const REAL_BDD *(*real_dd)(const EL_INFO *el_info, 1476 const QUAD *quad, int iq, void *apd); 1477 } LALt; 1478 MATENT_TYPE LALt_type; /* MATENT_REAL, _REAL_D or _REAL_DD */ 1479 bool LALt_pw_const; 1480 bool LALt_symmetric; 1481 int LALt_degree; 1482 1483 union { 1484 const REAL *(*real)(const EL_INFO *el_info, 1485 const QUAD *quad, int iq, void *apd); 1486 const REAL_D *(*real_d)(const EL_INFO *el_info, 1487 const QUAD *quad, int iq, void *apd); 1488 const REAL_DD *(*real_dd)(const EL_INFO *el_info, 1489 const QUAD *quad, int iq, void *apd); 1490 } Lb0; 1491 bool Lb0_pw_const; 1492 union { 1493 const REAL *(*real)(const EL_INFO *el_info, 1494 const QUAD *quad, int iq, void *apd); 1495 const REAL_D *(*real_d)(const EL_INFO *el_info, 1496 const QUAD *quad, int iq, void *apd); 1497 const REAL_DD *(*real_dd)(const EL_INFO *el_info, 1498 const QUAD *quad, int iq, void *apd); 1499 } Lb1; 1500 bool Lb1_pw_const; 1501 MATENT_TYPE Lb_type; /* MATENT_REAL, _REAL_D or _REAL_DD */ 1502 bool Lb0_Lb1_anti_symmetric; 1503 int Lb_degree; 1504 1505 union { 1506 REAL (*real)(const EL_INFO *el_info, 1507 const QUAD *quad, int iq, void *apd); 1508 const REAL *(*real_d)(const EL_INFO *el_info, 1509 const QUAD *quad, int iq, void *apd); 1510 const REAL_D *(*real_dd)(const EL_INFO *el_info, 1511 const QUAD *quad, int iq, void *apd); 1512 } c; 1513 bool c_pw_const; 1514 MATENT_TYPE c_type; /* MATENT_REAL, _REAL_D or _REAL_DD */ 1515 int c_degree; 1516 1517 BNDRY_FLAGS dirichlet_bndry; /* bndry-type bit-mask for 1518 * Dirichlet-boundary conditions 1519 * built into the matrix 1520 */ 1521 FLAGS fill_flag; 1522 void *user_data; /* application data, passed to init_element */ 1523}; 1524\end{lstlisting}\ev 1525 1526\begin{compatibility} 1527 Former versions of the \ALBERTA toolkit had special 1528 ``\code{DOWB\_OPERATOR\_INFO}'' and \code{DOF\_DOWB\_MATRIX}'' 1529 definitions to model block-matrix structures with $\DOW\times\DOW$ 1530 blocks, $1\times\DOW$ and $\DOW\times 1$ blocks and $1\times 1$ 1531 blocks (i.e. not-blocked). Because those structures included 1532 the scalar case as well, the ordinary scalar-only 1533 \code{OPERATOR\_INFO} and \code{DOF\_MATRIX} structures have been 1534 abandoned altogether, and the \code{\dots\_DOWB\_\dots} versions 1535 were renamed, dropping the bizarre \code{DOWB} component of their 1536 names. 1537\end{compatibility} 1538 1539Description of the \code{OPERATOR\_INFO} structure: 1540\begin{descr} 1541\hyperitem{OPERATOR_INFO:row_fe_space}{row\_fe\_space} pointer to a 1542 finite element space connected to the row DOFs of the resulting 1543 matrix. 1544 %% 1545\hyperitem{OPERATOR_INFO:col\_fe\_space}{col\_fe\_space} pointer to a 1546 finite element space connected to the column DOFs of the resulting 1547 matrix. 1548 %% 1549\hyperitem{OPERATOR_INFO:quad}{quad} vector with pointers to 1550 quadratures; \code{quad[0]} is used for the integration of the zero 1551 order term, \code{quad[1]} for the first order term(s), and 1552 \code{quad[2]} for the second order term. 1553 %% 1554 \idx{Per-element initializers!OPERATOR_INFO@{\code{OPERATOR\_INFO}}} 1555 \idx{Per-element initializers!BNDRY_OPERATOR_INFO@{\code{BNDRY\_OPERATOR\_INFO}}} 1556 \idx{init_element()@{\code{init\_element()}}!OPERATOR_INFO@{\code{OPERATOR\_INFO}}} 1557 \idx{init_element()@{\code{init\_element()}}!BNDRY_OPERATOR_INFO@{\code{BNDRY\_OPERATOR\_INFO}}} 1558\hyperitem{OPERATOR_INFO:init_element}{init\_element}\idx{init_element()@{\code{init\_element()}}} 1559 pointer to a function for doing an initialization step on each 1560 element; \code{init\_element} may be a \nil pointer; 1561 1562 if \code{init\_element} is not \nil, \code{init\_element(el\_info, 1563 quad, user\_data)} is the first statement executed on each element 1564 \code{el\_info->el} and may initialize data which is used by the 1565 functions \code{LALt()}, \code{Lb0()}, \code{Lb1()}, and/or 1566 \code{c()} (calculate the Jacobian of the barycentric coordinates in 1567 the 1st and 2nd order terms or the element volume for all order 1568 terms, e.g.); \code{quad} is a pointer to a vector of quadratures 1569 which is actually used for the integration of the various order 1570 terms and \code{user\_data} may hold a pointer to user data, filled 1571 by \code{init\_element()}, e.g.; the return value is of interest in 1572 the case of parametric meshes and should be \true if the element is 1573 a curved element and \false otherwise. 1574 %% 1575\hyperitem{OPERATOR_INFO:LALt}{LALt, LALt.real, LALt.real\_d, 1576 LALt.real\_dd}\idx{LALt()} is a pointer to a function for the 1577 evaluation of $\bar A$ at quadrature nodes on the element; 1578 \code{LALt} may be a \nil pointer, if no second order term has to be 1579 integrated. 1580 1581 if \code{LALt} is not \nil, \code{LALt(el\_info, quad, iq, 1582 user\_data)} returns a pointer to a matrix of size 1583 ${\code{N\_LAMBDA}}\times{\code{N\_LAMBDA}}$ storing the value of 1584 $\bar A$ at \code{quad->lambda[iq]}; \code{quad} is the quadrature 1585 for the second order term and \code{user\_data} is a pointer to user 1586 data and \code{EL\_INFO} the current element descriptor. 1587 1588 The element-type of the returned matrix is determined by 1589 \code{LALt\_type}, i.e. the actual return type is either 1590 \code{REAL\_BB} for \code{MATENT\_REAL}, \code{REAL\_BBD} for 1591 \code{MATENT\_REAL\_D} or \code{REAL\_BBDD} for 1592 \code{MATENT\_REAL\_DD}. Note that one of the \code{B}'s is missing 1593 in the structure definition above because the \code{LALt} is 1594 supposed to return the address of the first element of the matrix. 1595 %% 1596\hyperitem{OPERATOR_INFO:LALt_type}{LALt\_type} codes the 1597 block-matrix type, see \code{MATENT\_TYPE} on page 1598 \pageref{T:MATENT_TYPE}. 1599 %% 1600\hyperitem{OPERATOR_INFO:LALt_pw_const}{LALt\_pw\_const} should be 1601 \code{true} if $\bar A$ is piecewise constant on the mesh (constant 1602 matrix $A$ on a non--parametric mesh, e.g.); thus integration of the 1603 second order term can use pre--computed integrals of the basis 1604 functions on the standard element; otherwise integration is done by 1605 using quadrature on each element; this entry also influences the 1606 assembly on parametric meshes with \code{strategy>0}, see 1607 \secref{S:access_param_mesh}: \ALBERTA will assume a constant value 1608 of $\bar A$ for non--curved elements on a parametric mesh and 1609 optimize by only calling \code{LALt} once with \code{iq==0}; 1610 1611\hyperitem{OPERATOR_INFO:LALt_symmetric}{LALt\_symmetric} should be 1612 \code{true} if $\bar A$ is a symmetric matrix; if the finite element 1613 spaces for rows and columns are the same, only the diagonal and the 1614 upper part of the element matrix for the second order term have to 1615 be computed; elements of the lower part can then be set using the 1616 symmetry; otherwise the complete element matrix has to be 1617 calculated; 1618 1619\hyperitem{OPERATOR_INFO:LALt_degree}{LALt\_degree} If 1620 \code{LALt\_pw\_const == false}, the \code{LALt\_degree} gives a 1621 hint about which quadrature rule should be used to integrate the 1622 second order term. This has only an effect if 1623 \hyperlink{OPERATOR_INFO:quad}{\code{quad[2]} == \nil}. In that 1624 case, \ALBERTA takes \code{LALt\_degree} and the row- and column 1625 finite element spaces into account to select a suitable quadrature 1626 formula. 1627 1628\hyperitem{OPERATOR_INFO:Lb0}{Lb0, Lb0.real, Lb0.real\_d, 1629 Lb0.real\_dd} is a pointer to a function for the evaluation of 1630 $\bar b^0$, at quadrature nodes on the element; \code{Lb0} may be a 1631 \nil pointer, if this first order term has not to be integrated; 1632 1633 if \code{Lb0} is not \nil, \code{Lb0(el\_info, quad, iq, 1634 user\_data)} returns a pointer to a vector of length 1635 \code{N\_LAMBDA} storing the value of $\bar b^0$ at 1636 \code{quad->lambda[iq]}; \code{quad} is the quadrature for the first 1637 order term and \code{user\_data} is a pointer to user data; 1638 1639\hyperitem{OPERATOR_INFO:Lb0_pw_const}{Lb0\_pw\_const} should be 1640 \code{true} if $\bar b^0$ is piecewise constant on the mesh 1641 (constant vector $b$ on a non--parametric mesh, e.g.); thus 1642 integration of the first order term can use pre--computed integrals 1643 of the basis functions on the standard element; otherwise 1644 integration is done by using quadrature on each element; for 1645 parametric meshes the same remarks as for \code{LALt\_symmetric} 1646 above hold; 1647 1648\hyperitem{OPERATOR_INFO:Lb1}{Lb1, Lb1.real, Lb1.real\_d, 1649 Lb1.real\_dd} is a pointer to a function for the evaluation of 1650 $\bar b^1$, at quadrature nodes on the element; \code{Lb1} may be a 1651 \nil pointer, if this first order term has not to be integrated; 1652 1653 if \code{Lb1} is not \nil, \code{Lb1(el\_info, quad, iq, 1654 user\_data)} returns a pointer to a vector of length 1655 \code{N\_LAMBDA} storing the value of $\bar b^1$ at 1656 \code{quad->lambda[iq]}; \code{quad} is the quadrature for the first 1657 order term and \code{user\_data} is a pointer to user data; 1658 1659\hyperitem{OPERATOR_INFO:Lb1_pw_const}{Lb1\_pw\_const} should be 1660 \code{true} if $\bar b^1$ is piecewise constant on the mesh 1661 (constant vector $b$ on a non--parametric mesh, e.g.); thus 1662 integration of the first order term can use pre--computed integrals 1663 of the basis functions on the standard element; otherwise 1664 integration is done by using quadrature on each element; 1665 1666\hyperitem{OPERATOR_INFO:Lb_type}{Lb\_type} see \code{LALt\_type}. 1667 1668\hyperitem{OPERATOR_INFO:Lb0_Lb1_anti_symmetric}{Lb0\_Lb1\_anti\_symmetric} 1669 should be \code{true} if the contributions of the complete first 1670 order term to the \emph{local} element matrix are anti symmetric 1671 (only possible if both \code{Lb0} and \code{Lb1} are not \nil, $\bar 1672 b^0 = -\bar b^1$, e.g.); if the finite element spaces for rows and 1673 columns are the same then only the upper part of the element matrix 1674 for the first order term has to be computed; elements of the lower 1675 part can then be set using the anti symmetry; otherwise the complete 1676 element matrix has to be calculated; 1677 1678\hyperitem{OPERATOR_INFO:Lb_degree} See the explanations for 1679 \hyperlink{OPERATOR_INFO:LALt_degree}{\code{LALt\_degree}} above. 1680 1681\hyperitem{OPERATOR_INFO:c}{c, c.real, c.real\_d, c.real\_dd} is a 1682 pointer to a function for the evaluation of $\bar c$ at quadrature 1683 nodes on the element; \code{c} may be a \nil pointer, if no zero 1684 order term has to be integrated; 1685 1686 if \code{c} is not \nil, \code{c(el\_info, quad, iq, user\_data)} 1687 returns the value of the function $\bar c$ at 1688 \code{quad->lambda[iq]}; \code{quad} is the quadrature for the zero 1689 order term and \code{user\_data} is a pointer to user data; 1690 1691\hyperitem{OPERATOR_INFO:c_type}{c\_type} see \code{LALt\_type}. 1692 1693\hyperitem{OPERATOR_INFO:c_pw_const}{c\_pw\_const} should be 1694 \code{true} if the zero order term $\bar c$ is piecewise constant on 1695 the mesh (constant function $c$ on a non--parametric mesh, e.g.); 1696 thus integration of the zero order term can use pre--computed 1697 integrals of the basis functions on the standard element; otherwise 1698 integration is done by using quadrature on each element; the same 1699 remarks about parametric meshes as for the other \code{*\_pw\_const} 1700 entries hold; 1701 1702\hyperitem{OPERATOR_INFO:c_degree} See the explanations for 1703 \hyperlink{OPERATOR_INFO:LALt_degree}{\code{LALt\_degree}} above. 1704 1705\hyperitem{OPERATOR_INFO:dirichlet_bndry}{dirichlet\_bndry} A bit 1706 mask flagging those components of the boundary of the triangulation 1707 which are subject to Dirichlet boundary conditions. See 1708 \secref{S:boundary}. 1709 1710\hyperitem{OPERATOR_INFO:user_data}{user\_data} optional pointer to 1711 memory segment for ``user data'' used by \code{init\_element()}, 1712 \code{LALt()}, \code{Lb0()}, \code{Lb1()}, and/or \code{c()} and is 1713 the last argument to these functions. A better name would maybe 1714 ``application data'', at any rate this is the channel were an 1715 application program can communicate data -- like the determinant of 1716 the transformation to the reference element -- from 1717 \code{init\_element()} to the operator kernels \code{LALt} \& 1718 friends, without using global variables. \emph{The data behind this 1719 pointer must be persistent for the entire life time of the 1720 application program. Especially, \code{user\_data} must not point 1721 to the stack area of some sub-routine call.} 1722 1723\hyperitem{OPERATOR_INFO:fill_flag}{fill\_flag} the flag for the mesh 1724 traversal routine indicating which elements should be visited and 1725 which information should be present in the \code{EL\_INFO} structure 1726 for \code{init\_element()}, \code{LALt()}, \code{Lb0()}, 1727 \code{Lb1()}, and/or \code{c()} on the visited elements. 1728\end{descr} 1729 1730\hrulefill 1731 1732Sometimes it is necessary to add contributions of boundary integrals 1733to the system matrix. One example are ``Robin'' boundary conditions 1734(see \secref{S:robin_bound}), other important examples include 1735capillary boundary conditions in the context of free boundary 1736problems, or penalty terms to penalize tangential stresses. Another 1737context which requires integration over the boundaries of all mesh 1738elements is the implementation of discontinuous Galerkin (DG) methods. 1739To aid these tasks there is a \code{BNDRY\_OPERATOR\_INFO} structure, 1740which resembles in its layout the (bulk-) \code{OPERATOR\_INFO} 1741structure; it is defined as follows: 1742 1743\ddx{BNDRY_OPERATOR_INFO@{\code{BNDRY\_OPERATOR\_INFO}}} 1744\idx{assemblage tools!BNDRY_OPERATOR_INFO@{\code{BNDRY\_OPERATOR\_INFO}}} 1745%% 1746\newcommand{\BNDRYOPERATORINFO}{\hyperref[T:BNDRY_OPERATOR_INFO]{\code{BNDRY\_OPERATOR\_INFO}}\xspace} 1747%% 1748\bv\begin{lstlisting}[label=T:BNDRY_OPERATOR_INFO,name=BNDRY_OPERATOR_INFO] 1749typedef struct bndry_operator_info BNDRY_OPERATOR_INFO; 1750 1751struct bndry_operator_info 1752{ 1753 const FE_SPACE *row_fe_space; 1754 const FE_SPACE *col_fe_space; 1755 1756 const WALL_QUAD *quad[3]; 1757 1758 bool (*init_element)(const EL_INFO *el_info, int wall, 1759 const WALL_QUAD *quad[3], void *ud); 1760 1761 union { 1762 const REAL_B *(*real)(const EL_INFO *el_info, 1763 const QUAD *quad, int iq, void *apd); 1764 const REAL_BD *(*real_d)(const EL_INFO *el_info, 1765 const QUAD *quad, int iq, void *apd); 1766 const REAL_BDD *(*real_dd)(const EL_INFO *el_info, 1767 const QUAD *quad, int iq, void *apd); 1768 } LALt; 1769 MATENT_TYPE LALt_type; /* MATENT_REAL, _REAL_D or _REAL_DD */ 1770 bool LALt_pw_const; 1771 bool LALt_symmetric; 1772 int LALt_degree; 1773 1774 union { 1775 const REAL *(*real)(const EL_INFO *el_info, 1776 const QUAD *quad, int iq, void *apd); 1777 const REAL_D *(*real_d)(const EL_INFO *el_info, 1778 const QUAD *quad, int iq, void *apd); 1779 const REAL_DD *(*real_dd)(const EL_INFO *el_info, 1780 const QUAD *quad, int iq, void *apd); 1781 } Lb0; 1782 bool Lb0_pw_const; 1783 union { 1784 const REAL *(*real)(const EL_INFO *el_info, 1785 const QUAD *quad, int iq, void *apd); 1786 const REAL_D *(*real_d)(const EL_INFO *el_info, 1787 const QUAD *quad, int iq, void *apd); 1788 const REAL_DD *(*real_dd)(const EL_INFO *el_info, 1789 const QUAD *quad, int iq, void *apd); 1790 } Lb1; 1791 bool Lb1_pw_const; 1792 MATENT_TYPE Lb_type; /* MATENT_REAL, _REAL_D or _REAL_DD */ 1793 bool Lb0_Lb1_anti_symmetric; 1794 int Lb_degree; 1795 1796 union { 1797 REAL (*real)(const EL_INFO *el_info, 1798 const QUAD *quad, int iq, void *apd); 1799 const REAL *(*real_d)(const EL_INFO *el_info, 1800 const QUAD *quad, int iq, void *apd); 1801 const REAL_D *(*real_dd)(const EL_INFO *el_info, 1802 const QUAD *quad, int iq, void *apd); 1803 } c; 1804 bool c_pw_const; 1805 MATENT_TYPE c_type; /* MATENT_REAL, _REAL_D or _REAL_DD */ 1806 int c_degree; 1807 1808 /* boundary segment(s) we belong to; if 1809 * BNDRY_FLAGS_IS_INTERIOR(bndry_type), then the operator is invoked 1810 * on all interior faces, e.g. to implement a DG-method. 1811 */ 1812 BNDRY_FLAGS bndry_type; 1813 1814 bool discontinuous; /* assemble jumps w.r.t. the neighbour */ 1815 bool tangential; /* use tangential gradients */ 1816 1817 FLAGS fill_flag; 1818 void *user_data; 1819}; 1820\end{lstlisting}\ev 1821 1822Description: Because the general layout is the same as for the 1823bulk-\code{OPERATOR\_INFO} structure explained above we document only 1824the differences here. There are three additional components in the 1825structure: 1826 1827\begin{descr} 1828\hyperitem{BNDRY_OPERATOR_INFO:bndry_type}{bndry\_type} This is 1829 bit-mask and determines on which part of the boundary the operator 1830 should be invoked. See also \secref{S:boundary}. If 1831 \code{BNDRY\_FLAGS\_IS\_INTERIOR(bndry\_type)} evaluates to 1832 \code{true} (i.e. if bit $0$ is set in \code{bndry\_type}, then the 1833 operator is invoked on all walls of the triangulation, for instance 1834 to implement a DG-method. 1835 %% 1836\hyperitem{BNDRY_OPERATOR_INFO:discontinuous}{discontinuous} This is a 1837 boolean flag. If set to \code{true}, then the operator is treated as 1838 a DG-operator. This means, that it is invoked once for each wall of 1839 each element with the set of local basis functions on the neighbor 1840 element being used to define the column space (i.e. as 1841 ansatz-functions) and the set of local basis function on the current 1842 element defining the row-space (i.e. the test-functions). 1843 1844 One instance of \code{BNDRY\_OPERATOR\_INFO} can only be used to 1845 either implement a jump term or a term living on a single element. 1846 To have both, two instances have to be defined. To this aim 1847 \code{fill\_matrix\_info\_ext()} accepts multiple 1848 \code{BNDRY\_OPERATOR\_INFO} structures. The program-code 1849 \code{src/Common/ellipt-dg.c} in the \code{alberta-demo}-package 1850 implements a very simplistic DG-method as example: jumps over 1851 element boundaries are penalized by zero-order term. 1852 %% 1853\hyperitem{BNDRY_OPERATOR_INFO:tangential}{tangential} This is a 1854 boolean flag. If set to \code{true}, then only the tangential 1855 component of the gradients of the basis functions is used. 1856\end{descr} 1857 1858\hrulefill 1859 1860Information stored in \code{OPERATOR\_INFO} and 1861\code{BNDRY\_OPERATOR\_INFO} structures is used by the following 1862functions which return a pointer to a filled \code{EL\_MATRIX\_INFO} 1863structure; this structure can be used as an argument to the 1864\code{update\_matrix()} function which will then assemble the discrete 1865matrix corresponding to the operators defined by the 1866\code{OPERATOR\_INFO} and \code{BNDRY\_OPERATOR\_INFO} structures: 1867%% 1868 1869\fdx{fill_matrix_info()@{\code{fill\_matrix\_info()}}} 1870\idx{assemblage tools!fill_matrix_info()@{\code{fill\_matrix\_info()}}} 1871\fdx{fill_matrix_info_ext()@{\code{fill\_matrix\_info\_ext()}}} 1872\idx{assemblage tools!fill_matrix_info_)ext()@{\code{fill\_matrix\_info\_ext()}}} 1873\bv\begin{lstlisting} 1874const EL_MATRIX_INFO * 1875fill_matrix_info(const OPERATOR_INFO *operator_info, 1876 EL_MATRIX_INFO *matrix_info); 1877const EL_MATRIX_INFO * 1878fill_matrix_info_ext(EL_MATRIX_INFO *matrix_info, 1879 const OPERATOR_INFO *operator_info, 1880 const BNDRY_OPERATOR_INFO *bop_info, 1881 ... /* more bndry-ops */); 1882\end{lstlisting}\ev 1883Description: 1884\begin{descr} 1885\kitem{fill\_matrix\_info(op\_info, mat\_info)} 1886\kitem{fill\_matrix\_info\_ext(op\_info, mat\_info, bop\_info, \dots)} 1887 %% 1888 ~\hfill 1889 1890 Return a pointer to a filled \ELMATRIXINFO structure for the 1891 assemblage of the system matrix for the operator defined in 1892 \code{op\_info} and \code{bop\_info}. The difference between the two 1893 functions is that the \code{\dots\_ext()}-variant (``extended'') 1894 allows for additional arguments describing components of the 1895 differential operator which have to be assembled by boundary 1896 integrals. Multiple such boundary-operators can be passed to 1897 \code{fill\_matrix\_info\_ext()}, the final boundary operator must 1898 be followed by a \nil pointer. So 1899%% 1900 \bv\begin{lstlisting} 1901fill_matrix_info_ext(mat_info, operator_info, NULL); 1902\end{lstlisting}\ev 1903 is equivalent to 1904 \bv\begin{lstlisting} 1905fill_matrix_info(operator_info, mat_info); 1906\end{lstlisting}\ev 1907 There is the artificial restriction that at most 255 different 1908 \BNDRYOPERATORINFO structures may be passed. 1909 1910 If the argument \code{mat\_info} is a \nil pointer, a new structure 1911 \code{mat\_info} is allocated and filled; otherwise the structure 1912 \code{mat\_info} is filled; all members are newly assigned. 1913 1914 If the underlying finite element spaces form a direct sum, then this 1915 is taken care of automatically, and the return 1916 \ELMATRIXINFO structure will assemble block-matrices, 1917 where each block corresponds to the pairing of one component of the 1918 direct sum forming the ansatz-space and one component of the direct 1919 sum forming the space of test functions. See also 1920 \secref{S:chain_impl} and \secref{S:bfcts_chains} 1921 1922 The remaining part of this section is rather a description what 1923 happens ``back-stage'', when calling the 1924 \code{fill\_matrix\_info[\_ext]()} functions. The components of 1925 \ELMATRIXINFO are initialized like follows: 1926 1927 \begin{descr} 1928 \kitem{row\_fe\_space, col\_fe\_space} 1929 \code{op\_info->row\_fe\_space} and 1930 \code{op\_info->col\_fe\_space} are pointers to the finite element 1931 spaces (and by this to the basis functions and DOFs) connected to 1932 the row DOFs and the column DOFs of the matrix to be assembled; if 1933 both pointers are \nil pointers, an error message is given, and 1934 the program stops; if one of these pointers is \nil, rows and 1935 column DOFs are connected with the same finite element space (i.e. 1936 \code{op\_info->row\_fe\_space = op\_info->col\_fe\_space}, or 1937 \code{op\_info->col\_fe\_space = op\_info->row\_fe\_space} is 1938 used). 1939 1940 \kitem{krn\_blk\_type} Based on the matrix block-type of the zero, 1941 first and second order kernels 1942 \hyperlink{OPERATOR_INFO:LALt_type}{\code{oinfo->c\_type}}, 1943 \hyperlink{OPERATOR_INFO:LALt_type}{\code{oinfo->Lb\_type}} and 1944 \hyperlink{OPERATOR_INFO:LALt_type}{\code{oinfo->LALt\_type}} and 1945 on the dimension of the range of the row- and column finite 1946 element spaces \code{krn\_blk\_type} is set to either 1947 \code{MATENT\_REAL}, \code{MATENT\_REAL\_D} or 1948 \code{MATENT\_REAL\_DD} to reflect the block-matrix structure of 1949 the element matrix. 1950 1951 \kitem{dirichlet\_bndry} is just a copy of 1952 \hyperlink{OPERATOR_INFO:dirichlet_bndry}{\code{oinfo->dirichlet\_bndry}}, 1953 see also section \code{S:boundary}. 1954 1955 \kitem{factor} is initialized to \code{1.0}. Note that the structure 1956 returned carries the \code{const} qualifier; the clean way to 1957 obtain \ELMATRIXINFO structures with a modifiable \code{factor} 1958 component is to pass storage to \code{fill\_matrix\_info[\_ext]()} 1959 via the \code{matrix\_info} parameter. 1960 1961 \kitem{el\_matrix\_fct} The most important member in the structure, 1962 namely \code{mat\_info->el\_matrix\_fct}, is adjusted to some 1963 general routine for the integration of the element matrix for any 1964 set of local basis functions; \code{fill\_matrix\_info()} tries to 1965 use the fastest available function for the element integration for 1966 the operator defined in 1967 \hyperref[T:OPERATOR_INFO]{\code{op\_info}}, depending on 1968 \hyperlink{OPERATOR_INFO:LALt_pw_const}{\code{op\_info->LALt\_pw\_const}} 1969 and similar hints; 1970 1971 Denote by \code{row\_degree} and \code{col\_degree} the degree of 1972 the basis functions connected to the rows and columns. Internally, 1973 a three-element vector ``\code{quad}'' of quadratures rules is 1974 used for the element integration, if not specified by 1975 \code{op\_info->quad}. The quadratures are chosen according to the 1976 following rules: pre-computed integrals of basis functions should 1977 be evaluated exactly, and all terms calculated by quadrature on 1978 the elements should use the same quadrature formula (this is more 1979 efficient than to use different quadratures). To be more specific: 1980 1981 \begin{itemize} 1982 \item If the 2nd order term has to be integrated and 1983 \code{op\_info->quad[2]} is not \nil, \code{quad[2] = 1984 op\_info->quad[2]} is used, otherwise \code{quad[2]} is a 1985 quadrature which is exact of degree 1986 \code{row\_degree+col\_degree-2+oinfo->LALt\_degree}. If the 2nd 1987 order term is not integrated then \code{quad[2]} is set to \nil. 1988 1989 \item If the 1st order term has to be integrated and 1990 \code{op\_info->quad[1]} is not \nil, \code{quad[1] = 1991 op\_info->quad[1]} is used; otherwise: if 1992 \code{op\_info->Lb\_pw\_const} is zero and \code{quad[2]} is not 1993 \nil, \code{quad[1] = quad[2]} is used, otherwise \code{quad[1]} 1994 is a quadrature which is exact of degree 1995 \code{row\_degree+col\_degree-1+oinfo->Lb\_degree}. If the 1st 1996 order term is not integrated then \code{quad[1]} is set to \nil. 1997 1998 \item If the zero order term has to be integrated and 1999 \code{op\_info->quad[0]} is not \nil, \code{quad[0] = 2000 op\_info->quad[0]} is used; otherwise: if 2001 \code{op\_info->c\_pw\_const} is zero and \code{quad[2]} is not 2002 \nil, \code{quad[0] = quad[2]} is used, if \code{quad[2]} is 2003 \nil and \code{quad[1]} is not \nil, \code{quad[0] = quad[1]} is 2004 used, or if both quadratures are \nil, \code{quad[0]} is a 2005 quadrature which is exact of degree 2006 \code{row\_degree+col\_degree+oinfo->c\_degree}. If the zero 2007 order term is not integrated then \code{quad[0]} is set to \nil. 2008 \end{itemize} 2009 2010 \noindent 2011 \code{el\_matrix\_fct()} roughly works as follows: 2012 2013 \begin{itemize} 2014 \item If \code{op\_info->init\_element} is not \nil then a call of 2015 \code{op\_info->init\_element(el\_info, quad, 2016 op\_info->user\_data)} is the first statement of 2017 \code{mat\_info->el\_matrix\_fct()} on each element; 2018 \code{el\_info} is a pointer to the \code{EL\_INFO} structure of 2019 the actual element, \code{quad} is the quadrature vector 2020 described above (now giving information about the actually used 2021 quadratures), and the last argument is a pointer to the 2022 application-data pointer 2023 \hyperlink{OPERATOR_INFO:user_data}{\code{oinfo->user\_data}}. 2024 2025 \item If \code{op\_info->LALt} is not \nil, the 2nd order term is 2026 integrated using the quadrature \code{quad[2]}; if 2027 \code{op\_info->LALt\_pw\_const} is not zero, the integrals of 2028 the product of gradients of the basis functions on the standard 2029 simplex are initialized (using the quadrature \code{quad[2]} for 2030 the integration) and used for the computation on the elements; 2031 \code{op\_info->LALt()} is only called once with arguments 2032 \code{op\_info->LALt(el\_info, quad[2], 0, 2033 op\_info->user\_data)}, i.e. the matrix of the 2nd order term 2034 is evaluated only at the first quadrature node; otherwise the 2035 integrals are approximated by quadrature and 2036 \code{op\_info->LALt()} is called for each quadrature node of 2037 \code{quad[2]}; if \code{op\_info->LALt\_symmetric} is not zero, 2038 the symmetry of the element matrix is used, if the finite 2039 element spaces are the same and this term is not integrated by 2040 the same quadrature as the first order term. 2041 2042 \item If \code{op\_info->Lb0} is not \nil, this 1st order term is 2043 integrated using the quadrature \code{quad[1]}; if 2044 \code{op\_info->Lb0\_pw\_const} is not zero, the integrals of 2045 the product of basis functions with gradients of basis functions 2046 on the standard simplex are initialized (using the quadrature 2047 \code{quad[1]} for the integration) and used for the computation 2048 on the elements; \code{op\_info->Lb0()} is only called once with 2049 arguments \code{op\_info->Lb0(el\_info, quad[1], 0, 2050 op\_info->user\_data)}, i.e. the vector of this 1st order term 2051 is evaluated only at the first quadrature node; otherwise the 2052 integrals are approximated by quadrature and 2053 \code{op\_info->Lb0()} is called for each quadrature node of 2054 \code{quad[1]}; 2055 2056 \item If \code{op\_info->Lb1} is not \nil, this 1st order term is 2057 integrated also using the quadrature \code{quad[1]}; if 2058 \code{op\_info->Lb1\_pw\_const} is not zero, the integrals of 2059 the product of gradients of basis functions with basis functions 2060 on the standard simplex are initialized (using the quadrature 2061 \code{quad[1]} for the integration) and used for the computation 2062 on the elements; \code{op\_info->Lb1()} is only called once with 2063 arguments \code{op\_info->Lb1(el\_info, quad[1], 0, 2064 op\_info->user\_data)}, i.e. the vector of this 1st order term 2065 is evaluated only at the first quadrature node; otherwise the 2066 integrals are approximated by quadrature and 2067 \code{op\_info->Lb1()} is called for each quadrature node of 2068 \code{quad[1]}. 2069 2070 \item If both function pointers \code{op\_info->Lb0} and 2071 \code{op\_info->Lb1} are not \nil, the finite element spaces for 2072 rows and columns are the same and 2073 \code{Lb0\_Lb1\_anti\_symmetric} is non--zero, then the 2074 contributions of the 1st order term are computed using this anti 2075 symmetry property. 2076 2077 \item If \code{op\_info->c} is not \nil, the zero order term is 2078 integrated using the quadrature \code{quad[0]}; if 2079 \code{op\_info->c\_pw\_const} is not zero, the integrals of the 2080 product of basis functions on the standard simplex are 2081 initialized (using the quadrature \code{quad[0]} for the 2082 integration) and used for the computation on the elements; 2083 \code{op\_info->c()} is only called once with arguments 2084 \code{op\_info->c(el\_info, quad[0], 0, op\_info->user\_data)}, 2085 i.e. the zero order term is evaluated only at the first 2086 quadrature node; otherwise the integrals are approximated by 2087 quadrature and \code{op\_info->c()} is called for each 2088 quadrature node of \code{quad[0]}. 2089 2090 \item The functions \code{LALt()}, \code{Lb0()}, \code{Lb1()}, and 2091 \code{c()}, can be called in an arbitrary order on the elements, 2092 if not \nil (this depends on the type of integration, using 2093 pre--computed values, using same/different quadrature for the 2094 second, first, and/or zero order term, e.g.) but commonly used 2095 data for these functions is always initialized first by 2096 \code{op\_info->init\_element()}, if this function pointer is 2097 not \nil. 2098 2099 \item Using all information about the operator and quadrature, an 2100 ``optimal'' routine for the assemblage is chosen. Information 2101 for this routine is stored at \code{mat\_info} which includes 2102 the pointer to user data \code{op\_info->user\_data} (the last 2103 argument to \code{init\_element()}, \code{LALt()}, \code{Lb0()}, 2104 \code{Lb1()}, and/or \code{c()}). 2105 \end{itemize} 2106 2107 \kitem{neigh\_el\_mat\_fcts[]} See the documentation of the 2108 \hyperlink{BNDRY_OPERATOR_INFO:discontinuous}{discontinuous} 2109 component of the \BNDRYOPERATORINFO structure. 2110 2111 \kitem{fill\_flag} Finally, the flag for the mesh traversal used by 2112 the function \code{update\_matrix()} is set in 2113 \code{mat\_info->fill\_flag} to \code{op\_info->fill\_flag}; it 2114 indicates which elements should be visited and which information 2115 should be present in the \code{EL\_INFO} structure for 2116 \code{init\_element()}, \code{LALt()}, \code{Lb0/1()}, and/or 2117 \code{c()} on the visited elements. 2118 2119 If the boundary bit-mask \code{op\_info->dirichlet\_bndry} has 2120 bits set (see also \secref{S:boundary}), then the 2121 \code{FILL\_BOUND} flag is added to \code{mat\_info->fill\_flag}. 2122 \end{descr} 2123\end{descr} 2124 2125\begin{example}% 2126[Implementation of the differential operator {\boldmath$-\Delta$}]% 2127\label{Ex:LALt} 2128The following source fragment gives an example of the implementation 2129for the operator $-\Delta$ and the access to a \code{MATRIX\_INFO} 2130structure for the automatic assemblage of the system matrix for 2131this problem for any set of used basis functions. 2132 2133The source fragment shown here is part of the implementation for a 2134Poisson equation, which is the first model problem described in detail 2135in \secref{S:poisson-impl}. However, we will generalize the code given in 2136\secref{S:poisson-impl} to include the case of parametric meshes. 2137The assemblage of the discrete system including the load vector and 2138Dirichlet boundary values is spelled out in \secref{S:ellipt_build}. 2139 2140For the Poisson equation we only have to implement a constant second 2141order term. For passing information about the gradient of the 2142barycentric coordinates (at all quadrature points) from 2143\code{init\_element()} to the function 2144\code{LALt} we define the following structure 2145% 2146\bv\begin{lstlisting} 2147struct app_data 2148{ 2149 REAL_BD *Lambda; 2150 REAL *det; 2151}; 2152\end{lstlisting}\ev 2153% 2154The function \code{init\_element()} calculates the Jacobians $\Lambda$ 2155and determinants \code{det} of the barycentric coordinates and stores these 2156in the above defined 2157structure. In the case of a parametric mesh we fill \code{Lambda} with 2158the Jacobians and \code{det} with the determinants at all quadrature 2159points of \code{quad[2]}. For a non--parametric mesh we only fill the 2160zeroth entry of \code{Lambda} and \code{det}. If 2161\code{init\_element()} returns 2162\false, then \code{LALt()} is only called once for the current simplex with 2163\code{iq==0}, otherwise it is called for each quadrature point in 2164\code{quad[2]}. Note that we need a higher order quadrature than usual 2165to calculate the integral exactly for a curved parametric element. 2166% 2167\idx{init_element()@{\code{init\_element()}}!Example} 2168\bv\begin{lstlisting} 2169static bool init_element(const EL_INFO *el_info, const QUAD *quad[3], void *ud) 2170{ 2171 struct app_data *data = (struct app_data *)ud; 2172 PARAMETRIC *parametric = el_info->mesh->parametric; 2173 2174 if (parametric && parametric->init_element(el_info, parametric)) { 2175 parametric->grd_lambda(el_info, quad[2], 0, NULL, 2176 data->Lambda, NULL, data->det); 2177 return true; 2178 } else { 2179 data->det[0] = el_grd_lambda(el_info, data->Lambda[0]); 2180 return false; 2181 } 2182} 2183\end{lstlisting}\ev 2184%% 2185The function \code{LALt} now has to calculate the scaled matrix 2186product $|\det D F_S| \Lambda \Lambda^t$. Note that \code{LALt()} is 2187invoked only for the first quadrature point (\code{iq == 0}), if the 2188\code{OPERATOR\_INFO}-structure claims that the second-order kernel is 2189piece-wise constant and \code{parametric->init\_element()} returns 2190\code{false}, so using the index \code{iq} into the fields \code{det} 2191and \code{Lambda} does not access invalid data, and the assembling 2192linear systems remains relatively efficient, even in the context of 2193iso-parametric boundary approximation. 2194% 2195\idx{LALt()!parametric example} 2196\bv\begin{lstlisting} 2197const REAL_B *LALt(const EL_INFO *el_info, const QUAD *quad, 2198 int iq, void *ud) 2199{ 2200 struct app_data *data = (struct app_data *)ud; 2201 int i, j; 2202 static REAL_BB LALt; /* mind the "static" keyword! */ 2203 2204 for (i = 0; i < N_VERTICES(MESH_DIM); i++) { 2205 LALt[i][i] = SCP_DOW(data->Lambda[iq][i], data->Lambda[iq][i]); 2206 LALt[i][i] *= data->det[iq]; 2207 for (j = i+1; j < N_VERTICES(MESH_DIM); j++) { 2208 LALt[i][j] = SCP_DOW(data->Lambda[iq][i], data->Lambda[iq][j]); 2209 LALt[i][j] *= data->det[iq]; 2210 LALt[j][i] = LALt[i][j]; 2211 } 2212 } 2213 2214 return (const REAL_B *)LALt; 2215} 2216\end{lstlisting}\ev 2217 2218Before assembling the system matrix for the operator $-\Delta$, we 2219first have to initialize an \code{EL\_MATRIX\_INFO} structure. A 2220pointer to this \code{EL\_MATRIX\_INFO} structure is the second 2221argument to the function \code{update\_matrix()}, which finally 2222assembles the system matrix (compare~\secref{S:matrix_assemblage}). 2223 2224For the initialization we have to fill an \code{OPERATOR\_INFO} 2225structure collecting all information about the differential operator. 2226For $-\Delta$ we only need pointers to the functions 2227\code{init\_element()} and \code{LALt()} described above. The 2228differential operator is constant and symmetric, and information about 2229vertex coordinates is needed for computing the Jacobian of the 2230barycentric coordinates. Information about Dirichlet boundary values 2231should be assembled into the system matrix, hence the entry 2232\code{operator\_info->use\_get\_bound} is set \code{true}. 2233 2234The functions \code{init\_element()} and \code{LALt()} do not 2235depend on the finite element space which is used. This functions 2236can be used for any conforming finite element discretization for the 2237Poisson equation; all information needed about the actually used 2238finite elements is a pointer to this finite element space; we assume 2239that this pointer is stored in the variable \code{fe\_space}. 2240% 2241\bv\begin{lstlisting} 2242const EL_MATRIX_INFO *matrix_info = NULL; 2243static struct app_data app_data; /* Must be static! */ 2244OPERATOR_INFO o_info = { NULL, }; 2245 2246if(mesh->parametric) 2247 quad = get_quadrature(2, 2*fe_space->bas_fcts->degree + 2); 2248else 2249 quad = get_quadrature(2, 2*fe_space->bas_fcts->degree-2); 2250 2251app_data.Lambda = MEM_ALLOC(quad->n_points, REAL_BD); 2252app_data.det = MEM_ALLOC(quad->n_points, REAL); 2253 2254o_info.quad[2] = quad; 2255o_info.row_fe_space = o_info.col_fe_space = fe_space; 2256o_info.init_element = init_element; 2257o_info.LALt.real = LALt; 2258o_info.LALt_pw_const = true; /* pw const. assemblage is faster */ 2259o_info.LALt_symmetric = true; /* symmetric assemblage is faster */ 2260o_info.user_data = &app_data; /* application data */ 2261 2262/* Use, e.g., Dirichlet boundary conditions. */ 2263BNDRY_FLAGS_CPY(o_info.dirichlet_bndry, dirichlet_mask); 2264 2265o_info.fill_flag = CALL_LEAF_EL|FILL_COORDS; 2266 2267matrix_info = fill_matrix_info(&o_info, NULL); 2268\end{lstlisting}\ev 2269% 2270Full information about the operator is now available via the 2271\code{matrix\_info} structure and the system matrix \code{matrix} can then 2272easily be assembled for the selected finite element space by 2273\bv\begin{lstlisting} 2274 clear_dof_matrix(matrix); 2275 update_matrix(matrix, matrix_info, NoTranspose); 2276\end{lstlisting}\ev 2277\end{example} 2278 2279\subsection{Matrix assemblage for coupled second order problems}% 2280\label{S:matrix_assemblage_coupled}% 2281 2282The corresponding mechanism for coupled vector valued problems is very 2283similar, except for the additional indices necessary to describe 2284coupling. We start by stating the form of the element matrix with the 2285generalized first order term and different finite element spaces (see 2286also \ref{book:S:Dis2ndorder}): 2287 2288\begin{align*} 2289L_{S,\mu\nu}^{ij} &:= \int_{\Shat}\nablal 2290 \bar\psi^{i}(\lambda(\xhat)) \cdot \bar A^{\mu\nu}(\lambda(\xhat))\, 2291 \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat \\ 2292 &+ \int_{\Shat} \bar\psi^{i}(\lambda(\xhat))\; \bar 2293 b^{0,\mu\nu}(\lambda(\xhat)) \cdot \nablal \pbar^{j}(\lambda(\xhat)) 2294 \,d\xhat \\ 2295 &+ \int_{\Shat}^{} \nablal \bar\psi^{i}(\lambda(\xhat)) \cdot\; \bar 2296 b^{1,\mu\nu}(\lambda(\xhat))\, \pbar^{j}(\lambda(\xhat)) \,d\xhat\\ 2297 &+ \int_{\Shat} \bar c^{\mu\nu}(\lambda(\xhat))\, 2298 \bar\psi^{i}(\lambda(\xhat)) \, \pbar^{j}(\lambda(\xhat))\,d\xhat, 2299\end{align*} 2300with 2301\begin{align*} 2302 \bar A^{\mu\nu}(\lambda) &:= \left(\bar 2303 a_{kl}^{\mu\nu}(\lambda)\right)_{k,l = 0,\ldots,d} := |\det 2304 DF_S(\xhat(\lambda))| \, \Lambda(x(\lambda)) \, 2305 A^{\mu\nu}(x(\lambda)) \, \Lambda^t(x(\lambda)),\\ 2306 \bar b^{0,\mu\nu}(\lambda) & := \left(\bar 2307 b_{l}^{0,\mu\nu}(\lambda)\right)_{l = 0,\ldots,d} := |\det 2308 DF_S(\xhat(\lambda))| \, \Lambda(x(\lambda)) \, b^{0,\mu\nu}(x(\lambda)),\\ 2309 \bar b^{1,\mu\nu}(\lambda) & := \left(\bar 2310 b_{l}^{1,\mu\nu}(\lambda)\right)_{l = 0,\ldots,d} := |\det 2311 DF_S(\xhat(\lambda))| \, \Lambda(x(\lambda)) \, b^{1,\mu\nu}(x(\lambda)), 2312 \quad\text{and}\\ 2313 \bar c^{\mu\nu}(\lambda) & := |\det 2314 DF_S(\xhat(\lambda))| \, c^{\mu\nu}(x(\lambda)) 2315\end{align*} 2316for $\mu,\nu=1,\dots,n$, $n=\code{DIM\_OF\_WORLD}$. 2317 2318To store information about the coupled operator and finite element 2319spaces, we use the same \verb|OPERATOR_INFO| (see page 2320\pageref{T:OPERATOR_INFO}) structure as for the scalar problems, we 2321only have to adjust the respective \verb|MATENT_TYPE| structure 2322components to the correct block-matrix type. Also, the same 2323\verb|fill_matrix_info()| and \verb|add_element_matrix()| routines are 2324used for scalar and vector valued problems. 2325 2326\subsection{Data structures for storing pre-computed integrals of 2327 basis functions}\label{S:ass_info}% 2328\label{S:Q11_PSI_PHI} 2329\label{S:Q10_PSI_PHI} 2330\label{S:Q01_PSI_PHI} 2331\label{S:Q00_PSI_PHI} 2332 2333Assume a non--parametric triangulation and constant coefficient 2334functions $A$, $b$, and $c$. Since the Jacobian of the barycentric 2335coordinates is constant on $S$, the functions $\bar A_S$, $\bar b^0_S$, 2336$\bar b^1_S$, and $\bar c_S$ are constant on $S$ also. Now, looking at the 2337element matrix approximated by some quadrature $\hat Q$, we observe 2338\begin{equation}\label{e:psiphiformula} 2339\begin{split} 2340\hat Q\Big(\sum_{k,l = 0}^d (\bar a_{S,kl} \bar\psi^i_{,\lambda_k} \, 2341 \pbar^j_{,\lambda_l})\Big) 2342&= 2343\sum_{k,l = 0}^d \bar a_{S,kl} 2344\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j_{,\lambda_l}\Big), 2345\\ 2346\hat Q\Big(\sum_{l = 0}^d ( 2347\bar b^0_{S,l} \, \bar\psi^i \, \pbar^j_{,\lambda_l})\Big) 2348&= 2349\sum_{l = 0}^d \bar b^0_{S,l} \, 2350\hat Q\Big(\bar\psi^i \, \pbar^j_{,\lambda_l}\Big), 2351\\ 2352\hat Q\Big(\sum_{k = 0}^d ( 2353\bar b^1_{S,k} \, \bar\psi^i_{,\lambda_k} \, \pbar^j)\Big) 2354&= 2355\sum_{k = 0}^d \bar b^1_{S,k} \, 2356\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j)\Big), 2357\qquad\mbox{and}\\ 2358\hat Q\Big( 2359 (\bar c_S \, \bar\psi^i \,\pbar^j)\Big) 2360&= 2361\bar c_S\, \hat Q\Big(\bar\psi^i \,\pbar^j\Big). 2362\end{split} 2363\end{equation} 2364The values of the quadrature applied to the basis functions only 2365depend on the basis functions and the standard element but not on the 2366actual simplex $S$. All information about $S$ is given by $\bar A_S$, 2367$\bar b^0_S$, $\bar b^1_S$, and $\bar c_S$. Thus, these quadratures 2368have only to be calculated once, and can then be used on each element 2369during the assembling. 2370 2371For this we have to store for the basis functions 2372$\{\bar\psi_i\}_{i=1,\dots,n}$ and $\{\bar\vphi_j\}_{j=1,\dots,m}$ 2373the values 2374\begin{alignat*}{2} 2375\hat Q^{11}_{ij,kl} &:= 2376\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j_{,\lambda_l}\Big) 2377&\qquad &\mbox{for } 1\leq i \leq n,\; 1\leq j\leq m,\; 23780 \leq k,l \leq d, 2379\intertext{if $A \neq 0$,} 2380\hat Q^{01}_{ij,l} &:= 2381\hat Q\Big(\bar\psi^i \, \pbar^j_{,\lambda_l}\Big) 2382&&\mbox{for } 1\leq i \leq n,\; 1\leq j\leq m,\; 23830 \leq l \leq d, 2384\intertext{if $b^0 \neq 0$,} 2385\hat Q^{10}_{ij,k} &:= 2386\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j\Big) 2387&&\mbox{for } 1\leq i \leq n,\; 1\leq j\leq m,\; 23880 \leq k \leq d 2389\intertext{if $b^1 \neq 0$, and} 2390\hat Q^{00}_{ij} &:= 2391\hat Q\Big(\bar\psi^i \,\pbar^j\Big) 2392&&\mbox{for } 1\leq i \leq n,\; 1\leq j\leq m, 2393\end{alignat*} 2394if $c \neq 0$. Many of these values are zero, especially for the 2395first and second order terms (if $\bar\psi^i$ and $\pbar^j$ are the 2396linear nodal basis functions $\hat Q^{11}_{ij,kl} = \delta_{ij}\delta_{kl}$). 2397Thus, we use special data structures 2398for a sparse storage of the non zero values for these terms. 2399These are described now. 2400 2401In order to ``define'' zero entries we use 2402\bv\begin{lstlisting} 2403static const REAL TOO_SMALL = 10.0 * REAL_EPSILON; 2404\end{lstlisting}\ev 2405and all computed values \code{val} with 2406$|\code{val}|\leq\code{TOO\_SMALL}$ are treated as zeros. As we are 2407considering here integrals over the standard simplex, non-zero 2408integrals are usually of order one, such that the above constant is of 2409the order of roundoff errors for double precision. 2410 2411The following data structure is used for storing values $\hat Q^{11}$ 2412for two sets of basis functions integrated with a given quadrature. 2413Note that in the context of ``chained'' basis-functions (see 2414\secref{S:bfcts_chains} the cache-structure nevertheless hold data for 2415only a single component of such a multi-component chain. 2416 2417\ddx{Q11_PSI_PHI@{\code{Q11\_PSI\_PHI}}}% 2418\idx{assemblage tools!Q11_PSI_PHI@{\code{Q11\_PSI\_PHI}}}% 2419\ddx{Q11_PSI_PHI_CACHE@{\code{Q11\_PSI\_PHI\_CACHE}}}% 2420\idx{assemblage tools!Q11_PSI_PHI_CHACHE@{\code{Q11\_PSI\_PHI\_CACHE}}}% 2421\bv\begin{lstlisting}[label=T:Q11_PSI_PHI] 2422typedef struct q11_psi_phi Q11_PSI_PHI; 2423 2424struct q11_psi_phi 2425{ 2426 const BAS_FCTS *psi; 2427 const BAS_FCTS *phi; 2428 const QUAD *quad; 2429 2430 const Q11_PSI_PHI_CACHE *cache; 2431 2432 INIT_ELEMENT_DECL; 2433}; 2434 2435typedef struct q11_psi_phi_cache 2436{ 2437 int n_psi; 2438 int n_phi; 2439 2440 const int *const*n_entries; 2441 const REAL *const*const*values; 2442 const int *const*const*k; 2443 const int *const*const*l; 2444} Q11_PSI_PHI_CACHE; 2445\end{lstlisting}\ev 2446Description: 2447\begin{descr} 2448 %% 2449 \kitem{struct q11\_psi\_phi}:\hfill 2450 \begin{descr} 2451 %% 2452 \kitem{psi} Pointer to the first set of basis functions. 2453 %% 2454 \kitem{phi} Pointer to the second set of basis functions. 2455 %% 2456 \kitem{quad} Pointer to the quadrature which is used for the integration. 2457 %% 2458 \kitem{cache} Pointer to the actual data in the cache. 2459 %% 2460 \kitem{INIT\_ELEMENT\_DECL} Optional per-element initializer. This 2461 entry is initialized when calling \code{get\_q11\_psi\_phi()} if 2462 either the underlying basis functions or the underlying quadrature 2463 rule has per-element initializers. See \secref{S:init_element}. 2464 %% 2465 \end{descr} 2466 \kitem{struct q11\_psi\_phi\_cache}:\hfill\ 2467 \begin{descr} 2468 \kitem{n\_psi} Dimension of the local space of test-functions (row 2469 space), equals \code{Q11\_PSI\_PHI.psi->n\_bas\_fcts}. 2470 %% 2471 \kitem{n\_phi} Dimension of the local space of ansatz-functions 2472 (column space), equals \code{Q11\_PSI\_PHI.phi->n\_bas\_fcts}. 2473 %% 2474 \kitem{n\_entries} matrix of size $\code{n\_psi}\times\code{n\_phi}$ 2475 storing the count of non zero integrals; 2476 2477 \code{n\_entries[i][j]} is the count of non zero values of $\hat 2478 Q^{11}_{\code{ij,kl}}$ ($0\leq\code{k,l}\leq d$) for the pair 2479 $(\code{psi[i]},\code{phi[j]})$, 2480 $\code{0}\leq\code{i}<\code{n\_psi}$, 2481 $\code{0}\leq\code{j}<\code{n\_phi}$. 2482 % 2483 \kitem{values} tensor storing the non zero integrals; 2484 2485 \code{values[i][j]} is a vector of length \code{n\_entries[i][j]} 2486 storing the non zero values for the pair 2487 $(\code{psi[i]},\code{phi[j]})$. 2488 %% 2489 \kitem{k, l} tensor storing the indices $k$, $l$ of the non zero 2490 integrals; 2491 2492 \code{k[i][j]} and \code{l[i][j]} are vectors of length 2493 \code{n\_entries[i][j]} storing at \code{k[i][j][r]} and 2494 \code{l[i][j][r]} the indices \code{k} and \code{l} of the value 2495 stored at \code{values[i][j][r]}, i.e. 2496 \end{descr} 2497\end{descr} 2498The following formulas summarize the relationship between the cache 2499data-structure and the formulas \eqref{e:psiphiformula} at the 2500beginning of this section: 2501\[ 2502\code{values[i][j][r]} = 2503\hat Q^{11}_{\code{ij},\code{k[i][j][r],l[i][j][r]}} 2504= \hat Q\Big(\bar\psi^{\code{i}}_{,\lambda_{\code{k[i][j][r]}}} \, 2505\pbar^{\code{j}}_{,\lambda_{\code{l[i][j][r]}}}\Big), 2506\] 2507for $\code{0}\leq\code{r}<\code{n\_entries[i][j]}$. 2508Using these pre--computed values we have for all elements $S$ 2509\[ 2510\sum_{k,l = 0}^d \bar a_{S,kl} 2511\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j_{,\lambda_l}\Big) 2512= 2513\sum_{\code{r}=\code{0}}^{\code{n\_entries[i][j]-1}} 2514\bar a_{S,\code{k[i][j][r],l[i][j][r]}}\,\code{*} 2515 {\code{values[i][j][r]}}. 2516\] 2517The following function initializes a \code{Q11\_PSI\_PHI} structure: 2518\fdx{get_q11_psi_phi()@{\code{get\_q11\_psi\_phi())}}}% 2519\idx{assemblage tools!get_q11_psi_phi()@{\code{get\_q11\_psi\_phi()}}}% 2520\bv\begin{lstlisting} 2521const Q11_PSI_PHI *get_q11_psi_phi(const BAS_FCTS *psi, const BAS_FCTS *phi, 2522 const QUAD *quad); 2523\end{lstlisting}\ev 2524Description: 2525\begin{descr} 2526 \kitem{get\_q11\_psi\_phi(psi, phi, quad)} returns a pointer to a filled 2527 \code{Q11\_PSI\_PHI} structure. 2528 2529 \code{psi} is a pointer to the first set of basis functions, 2530 \code{phi} is a pointer to the second set of basis functions; 2531 if both are \nil pointers, nothing is done and the return value 2532 is \nil; if one of the pointers is a \nil pointer, the structure 2533 is initialized using the same set of basis functions for the 2534 first and second set, i.e. \code{phi = psi} or \code{psi = phi} 2535 is used. 2536 2537 \code{quad} is a pointer to a quadrature for the approximation 2538 of the integrals; if \code{quad} is \nil, then a quadrature which 2539 is exact of degree \code{psi->degree+phi->degree-2} is used. 2540 2541 All used \code{Q11\_PSI\_PHI} structures are stored in a linked 2542 list and are identified uniquely by the members \code{psi}, 2543 \code{phi}, and \code{quad}, and for such a combination only 2544 one \code{Q11\_PSI\_PHI} structure is created during runtime; 2545 2546 First, \code{get\_q11\_psi\_phi()} 2547 looks for a matching structure in the linked list; if such a 2548 structure is found a pointer to this structure is returned; 2549 the values are not computed a second time. 2550 Otherwise a new structure is generated, linked to the list, and 2551 the values are computed using the quadrature \code{quad}; all 2552 values \code{val} with $|\code{val}|\leq\code{TOO\_SMALL}$ are 2553 treated as zeros. 2554\end{descr} 2555 2556\begin{example} 2557The following example shows how to use these pre--computed values 2558for the integration of the 2nd order term 2559\[ 2560\int_{\Shat} 2561 \nablal \bar\psi^{i}(\lambda(\xhat)) \cdot \bar A(\lambda(\xhat))\, 2562 \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat 2563\] 2564for all $i,j$. We only show the body of a function for the integration 2565and assume that \code{LALt\_fct} returns a matrix storing $\bar A$ 2566(compare the member \code{LALt} in the \code{OPERATOR\_INFO} structure): 2567\bv\begin{lstlisting} 2568 static Q11_PSI_PHI_CACHE *q11_psi_phi; 2569 2570 if (!q11_psi_phi) { 2571 q11_psi_phi = get_q11_psi_phi(psi, phi, quad[2])->cache; 2572 } 2573 2574 LALt = LALt_fct(el_info, quad, 0, user_data); 2575 n_entries = q11_psi_phi->n_entries; 2576 2577 for (i = 0; i < q11_psi_phi->n_psi; i++) 2578 { 2579 for (j = 0; j < q11_psi_phi->n_phi; j++) 2580 { 2581 k = q11_psi_phi->k[i][j]; 2582 l = q11_psi_phi->l[i][j]; 2583 values = q11_psi_phi->values[i][j]; 2584 for (val = m = 0; m < n_entries[i][j]; m++) 2585 val += values[m]*LALt[k[m]][l[m]]; 2586 mat[i][j] += val; 2587 } 2588 } 2589\end{lstlisting}\ev 2590\end{example} 2591 2592The values $\hat Q^{01}$ for the set of basis functions \code{psi} 2593and \code{phi} are stored in 2594\ddx{Q01_PSI_PHI@{\code{Q01\_PSI\_PHI}}}% 2595\idx{assemblage tools!Q01_PSI_PHI@{\code{Q01\_PSI\_PHI}}}% 2596\ddx{Q01_PSI_PHI_CACHE@{\code{Q01\_PSI\_PHI\_CACHE}}}% 2597\idx{assemblage tools!Q01_PSI_PHI_CACHE@{\code{Q01\_PSI\_PHI\_CACHE}}}% 2598\bv\begin{lstlisting}[label=T:Q01_PSI_PHI] 2599typedef struct q01_psi_phi Q01_PSI_PHI; 2600 2601typedef struct q01_psi_phi_cache 2602{ 2603 int n_psi; 2604 int n_phi; 2605 2606 const int *const*n_entries; 2607 const REAL *const*const*values; 2608 const int *const*const*l; 2609} Q01_PSI_PHI_CACHE; 2610 2611struct q01_psi_phi 2612{ 2613 const BAS_FCTS *psi; 2614 const BAS_FCTS *phi; 2615 const QUAD *quad; 2616 2617 const Q01_PSI_PHI_CACHE *cache; 2618 2619 INIT_ELEMENT_DECL; 2620}; 2621\end{lstlisting}\ev 2622Description: 2623\begin{descr} 2624 \kitem{struct q01\_psi\_phi}:\hfill 2625 \begin{descr} 2626 %% 2627 \kitem{psi} pointer to the first set of basis functions. 2628 %% 2629 \kitem{phi} pointer to the second set of basis functions. 2630 %% 2631 \kitem{quad} pointer to the quadrature which is used for the integration. 2632 %% 2633 \kitem{cache} Pointer to the actual data in the cache. 2634 %% 2635 \kitem{INIT\_ELEMENT\_DECL} Optional per-element initializer. This 2636 entry is initialized when calling \code{get\_q11\_psi\_phi()} if 2637 either the underlying basis functions or the underlying quadrature 2638 rule has per-element initializers. See \secref{S:init_element}. 2639 %% 2640 \end{descr} 2641 \kitem{struct q01\_psi\_phi\_cache}:\hfill 2642 \begin{descr} 2643 \kitem{n\_psi} Dimension of the local space of test-functions (row 2644 space), equals \code{Q11\_PSI\_PHI.psi->n\_bas\_fcts}. 2645 %% 2646 \kitem{n\_phi} Dimension of the local space of ansatz-functions 2647 (column space), equals \code{Q11\_PSI\_PHI.phi->n\_bas\_fcts}. 2648 %% 2649 \kitem{n\_entries} matrix of size 2650 $\code{psi->n\_bas\_fcts}\times\code{phi->n\_bas\_fcts}$ 2651 storing the count of non zero integrals; 2652 2653 \code{n\_entries[i][j]} is the count of non zero values 2654 of $\hat Q^{01}_{\code{ij,l}}$ ($0\leq\code{l}\leq d$) 2655 for the pair $(\code{psi[i]},\code{phi[j]})$, 2656 $\code{0}\leq\code{i}<\code{n\_psi}$, 2657 $\code{0}\leq\code{j}<\code{n\_phi}$. 2658 %% 2659 \kitem{values} tensor storing the non zero integrals; 2660 2661 \code{values[i][j]} is a vector of length \code{n\_entries[i][j]} 2662 storing the non zero values for the pair 2663 $(\code{psi[i]},\code{phi[j]})$. 2664 %% 2665 \kitem{l} tensor storing the indices $l$ of the non zero integrals; 2666 2667 \code{l[i][j]} is a vector of length \code{n\_entries[i][j]} 2668 storing at \code{l[i][j][r]} the index \code{l} of the 2669 value stored at \code{values[i][j][r]}, i.e. 2670 %% 2671 \end{descr} 2672\end{descr} 2673The following formulas summarize the relationship between the cache 2674data-structure and the formulas \eqref{e:psiphiformula} at the 2675beginning of this section: 2676\[ 2677\code{values[i][j][r]} = 2678\hat Q^{01}_{\code{ij},\code{l[i][j][r]}} 2679= \hat Q\Big(\bar\psi^{\code{i}} \, 2680\pbar^{\code{j}}_{,\lambda_{\code{l[i][j][r]}}}\Big), 2681\] 2682for $\code{0}\leq\code{r}<\code{n\_entries[i][j]}$. 2683Using these pre--computed values we have for all elements $S$ 2684\[ 2685\sum_{l = 0}^d \bar b^0_{S,l} 2686\hat Q\Big(\bar\psi^i \, \pbar^j_{,\lambda_l}\Big) 2687= 2688\sum_{\code{r}=\code{0}}^{\code{n\_entries[i][j]-1}} 2689\bar b^0_{S,\code{l[i][j][r]}}\,\code{*} {\code{values[i][j][r]}}. 2690\] 2691The following function initializes a \code{Q01\_PSI\_PHI} structure: 2692\fdx{get_q01_psi_phi()@{\code{get\_q01\_psi\_phi())}}}% 2693\idx{assemblage tools!get_q01_psi_phi()@{\code{get\_q01\_psi\_phi()}}}% 2694\bv\begin{lstlisting} 2695const Q01_PSI_PHI *get_q01_psi_phi(const BAS_FCTS *psi, const BAS_FCTS *phi, 2696 const QUAD *quad); 2697\end{lstlisting}\ev 2698Description: 2699\begin{descr} 2700\kitem{get\_q01\_psi\_phi(psi, phi, quad)} returns a pointer to a filled 2701 \code{Q01\_PSI\_PHI} structure. 2702 2703 \code{psi} is a pointer to the first set of basis functions 2704 \code{phi} is a pointer to the second set of basis functions; 2705 if both are \nil pointers, nothing is done and the return value 2706 is \nil; if one of the pointers is a \nil pointer, the structure 2707 is initialized using the same set of basis functions for the 2708 first and second set, i.e. \code{phi = psi} or \code{psi = phi} 2709 is used. 2710 2711 \code{quad} is a pointer to a quadrature for the approximation 2712 of the integrals; is \code{quad} is \nil, a quadrature which 2713 is exact of degree \code{psi->degree+phi->degree-1} is used. 2714 2715 All used \code{Q01\_PSI\_PHI} structures are stored in a linked 2716 list and are identified uniquely by the members \code{psi}, 2717 \code{phi}, and \code{quad}, and for such a combination only 2718 one \code{Q01\_PSI\_PHI} structure is created during runtime; 2719 2720 First, \code{get\_q01\_psi\_phi()} 2721 looks for a matching structure in the linked list; if such a 2722 structure is found a pointer to this structure is returned; 2723 the values are not computed a second time. 2724 Otherwise a new structure is generated, linked to the list, and 2725 the values are computed using the quadrature \code{quad}; all 2726 values \code{val} with $|\code{val}|\leq\code{TOO\_SMALL}$ are 2727 treated as zeros. 2728\end{descr} 2729The values $\hat Q^{10}$ for the set of basis functions \code{psi} 2730and \code{phi} are stored in 2731\ddx{Q10_PSI_PHI@{\code{Q10\_PSI\_PHI}}}% 2732\idx{assemblage tools!Q10_PSI_PHI@{\code{Q10\_PSI\_PHI}}}% 2733\ddx{Q10_PSI_PHI_CACHE@{\code{Q10\_PSI\_PHI\_CACHE}}}% 2734\idx{assemblage tools!Q10_PSI_PHI_CACHE@{\code{Q10\_PSI\_PHI\_CACHE}}}% 2735\bv\begin{lstlisting}[label=T:Q10_PSI_PHI] 2736typedef struct q10_psi_phi Q10_PSI_PHI; 2737 2738typedef struct q10_psi_phi_cache 2739{ 2740 int n_psi; 2741 int n_phi; 2742 2743 const int *const*n_entries; 2744 const REAL *const*const*values; 2745 const int *const*const*k; 2746} Q10_PSI_PHI_CACHE; 2747 2748struct q10_psi_phi 2749{ 2750 const BAS_FCTS *psi; 2751 const BAS_FCTS *phi; 2752 const QUAD *quad; 2753 2754 const Q10_PSI_PHI_CACHE *cache; 2755 2756 INIT_ELEMENT_DECL; 2757}; 2758\end{lstlisting}\ev 2759Description: 2760\begin{descr} 2761 \kitem{struct q10\_psi\_phi}:\hfill 2762 \begin{descr} 2763 \kitem{psi} pointer to the first set of basis functions. 2764 %% 2765 \kitem{phi} pointer to the second set of basis functions. 2766 %% 2767 \kitem{quad} pointer to the quadrature which is used for the integration. 2768 %% 2769 %% 2770 \kitem{cache} Pointer to the actual data in the cache. 2771 %% 2772 \kitem{INIT\_ELEMENT\_DECL} Optional per-element initializer. This 2773 entry is initialized when calling \code{get\_q11\_psi\_phi()} if 2774 either the underlying basis functions or the underlying quadrature 2775 rule has per-element initializers. See \secref{S:init_element}. 2776 %% 2777 \end{descr} 2778 \kitem{struct q10\_psi\_phi\_cache}:\hfill 2779 \begin{descr} 2780 \kitem{n\_psi} Dimension of the local space of test-functions (row 2781 space), equals \code{Q11\_PSI\_PHI.psi->n\_bas\_fcts}. 2782 %% 2783 \kitem{n\_phi} Dimension of the local space of ansatz-functions 2784 (column space), equals \code{Q11\_PSI\_PHI.phi->n\_bas\_fcts}. 2785 %% 2786 \kitem{n\_entries} matrix of size 2787 $\code{psi->n\_bas\_fcts}\times\code{phi->n\_bas\_fcts}$ 2788 storing the count of non zero integrals; 2789 2790 \code{n\_entries[i][j]} is the count of non zero values 2791 of $\hat Q^{10}_{\code{ij,k}}$ ($0\leq\code{k}\leq d$) 2792 for the pair $(\code{psi[i]},\code{phi[j]})$, 2793 $\code{0}\leq\code{i}<\code{n\_psi}$, 2794 $\code{0}\leq\code{j}<\code{n\_phi}$. 2795 %% 2796 \kitem{values} tensor storing the non zero integrals; 2797 2798 \code{values[i][j]} is a vector of length \code{n\_entries[i][j]} 2799 storing the non zero values for the pair 2800 $(\code{psi[i]},\code{phi[j]})$. 2801 %% 2802 \kitem{k} tensor storing the indices $k$ of the non zero integrals; 2803 2804 \code{k[i][j]} is a vector of length \code{n\_entries[i][j]} 2805 storing at \code{k[i][j][r]} the index \code{k} of the 2806 value stored at \code{values[i][j][r]}, i.e. 2807 %% 2808 \end{descr} 2809\end{descr} 2810The following formulas summarize the relationship between the cache 2811data-structure and the formulas \eqref{e:psiphiformula} at the 2812beginning of this section: 2813\[ 2814\code{values[i][j][r]} = 2815\hat Q^{10}_{\code{ij},\code{k[i][j][r]}} 2816= \hat Q\Big(\bar\psi^{\code{i}}_{,\lambda_{\code{k[i][j][r]}}} \, 2817\pbar^{\code{j}}\Big), 2818\] 2819for $\code{0}\leq\code{r}<\code{n\_entries[i][j]}$. 2820Using these pre--computed values we have for all elements $S$ 2821\[ 2822\sum_{k = 0}^d \bar b^1_{S,k} 2823\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j\Big) 2824= 2825\sum_{\code{r}=\code{0}}^{\code{n\_entries[i][j]-1}} 2826\bar b^1_{S,\code{k[i][j][r]}}\,\code{*} {\code{values[i][j][r]}}. 2827\] 2828The following function initializes a \code{Q10\_PSI\_PHI} structure: 2829\fdx{get_q10_psi_phi()@{\code{get\_q10\_psi\_phi())}}}% 2830\idx{assemblage tools!get_q10_psi_phi()@{\code{get\_q10\_psi\_phi()}}}% 2831\bv\begin{lstlisting} 2832const Q10_PSI_PHI *get_q10_psi_phi(const BAS_FCTS *psi, const BAS_FCTS *phi, 2833 const QUAD *quad); 2834\end{lstlisting}\ev 2835Description: 2836\begin{descr} 2837\kitem{get\_q10\_psi\_phi(psi, phi, quad)} returns a pointer to a filled 2838 \code{Q10\_PSI\_PHI} structure. 2839 2840 \code{psi} is a pointer to the first set of basis functions 2841 \code{phi} is a pointer to the second set of basis functions; 2842 if both are \nil pointers, nothing is done and the return value 2843 is \nil; if one of the pointers is a \nil pointer, the structure 2844 is initialized using the same set of basis functions for the 2845 first and second set, i.e. \code{phi = psi} or \code{psi = phi} 2846 is used. 2847 2848 \code{quad} is a pointer to a quadrature for the approximation 2849 of the integrals; is \code{quad} is \nil, a quadrature which 2850 is exact of degree \code{psi->degree+phi->degree-1} is used. 2851 2852 All used \code{Q10\_PSI\_PHI} structures are stored in a linked 2853 list and are identified uniquely by the members \code{psi}, 2854 \code{phi}, and \code{quad}, and for such a combination only 2855 one \code{Q10\_PSI\_PHI} structure is created during runtime; 2856 2857 First, \code{get\_q10\_psi\_phi()} 2858 looks for a matching structure in the linked list; if such a 2859 structure is found a pointer to this structure is returned; 2860 the values are not computed a second time. 2861 Otherwise a new structure is generated, linked to the list, and 2862 the values are computed using the quadrature \code{quad}; all 2863 values \code{val} with $|\code{val}|\leq\code{TOO\_SMALL}$ are 2864 treated as zeros. 2865\end{descr} 2866Finally, the values $\hat Q^{00}$ for the set of basis functions \code{psi} 2867and \code{phi} are stored in 2868\ddx{Q00_PSI_PHI@{\code{Q00\_PSI\_PHI}}}% 2869\idx{assemblage tools!Q00_PSI_PHI@{\code{Q00\_PSI\_PHI}}}% 2870\ddx{Q00_PSI_PHI_CACHE@{\code{Q00\_PSI\_PHI\_CACHE}}}% 2871\idx{assemblage tools!Q00_PSI_PHI_CACHE@{\code{Q00\_PSI\_PHI\_CACHE}}}% 2872\bv\begin{lstlisting}[label=T:Q00_PSI_PHI] 2873typedef struct q00_psi_phi Q00_PSI_PHI; 2874 2875typedef struct q00_psi_phi_cache 2876{ 2877 int n_psi; 2878 int n_phi; 2879 2880 const REAL *const*values; 2881} Q00_PSI_PHI_CACHE; 2882 2883struct q00_psi_phi 2884{ 2885 const BAS_FCTS *psi; 2886 const BAS_FCTS *phi; 2887 const QUAD *quad; 2888 2889 const Q00_PSI_PHI_CACHE *cache; 2890 2891 INIT_ELEMENT_DECL; 2892}; 2893 2894\end{lstlisting}\ev 2895Description: 2896\begin{descr} 2897 \kitem{struct q00\_psi\_phi}:\hfill 2898 \begin{descr} 2899 \kitem{psi} pointer to the first set of basis functions. 2900 %% 2901 \kitem{phi} pointer to the second set of basis functions. 2902 %% 2903 \kitem{quad} pointer to the quadrature which is used for the integration. 2904 %% 2905 %% 2906 \kitem{cache} Pointer to the actual data in the cache. 2907 %% 2908 \kitem{INIT\_ELEMENT\_DECL} Optional per-element initializer. This 2909 entry is initialized when calling \code{get\_q11\_psi\_phi()} if 2910 either the underlying basis functions or the underlying quadrature 2911 rule has per-element initializers. See \secref{S:init_element}. 2912 %% 2913 \end{descr} 2914 \kitem{struct q00\_psi\_phi\_cache}:\hfill 2915 \begin{descr} 2916 \kitem{n\_psi} Dimension of the local space of test-functions (row 2917 space), equals \code{Q11\_PSI\_PHI.psi->n\_bas\_fcts}. 2918 %% 2919 \kitem{n\_phi} Dimension of the local space of ansatz-functions 2920 (column space), equals \code{Q11\_PSI\_PHI.phi->n\_bas\_fcts}. 2921 %% 2922 \kitem{values} matrix storing the integrals. 2923 \end{descr} 2924\end{descr} 2925The following formulas summarize the relationship between the cache 2926data-structure and the formulas \eqref{e:psiphiformula} at the 2927beginning of this section: 2928\[ 2929\code{values[i][j]} = 2930\hat Q^{00}_{\code{ij}} 2931= \hat Q\Big(\bar\psi^{\code{i}} \, \pbar^{\code{j}}\Big), 2932\] 2933 for the pair $(\code{psi[i]},\code{phi[j]})$, 2934 $\code{0}\leq\code{i}<\code{psi->n\_bas\_fcts}$, 2935 $\code{0}\leq\code{j}<\code{phi->n\_bas\_fcts}$. 2936The following function initializes a \code{Q00\_PSI\_PHI} structure: 2937\fdx{get_q00_psi_phi()@{\code{get\_q00\_psi\_phi())}}}% 2938\idx{assemblage tools!get_q00_psi_phi()@{\code{get\_q00\_psi\_phi()}}}% 2939\bv\begin{lstlisting} 2940const Q00_PSI_PHI *get_q00_psi_phi(const BAS_FCTS *psi, const BAS_FCTS *phi, 2941 const QUAD *quad); 2942\end{lstlisting}\ev 2943Description: 2944\begin{descr} 2945\kitem{get\_q00\_psi\_phi(psi, phi, quad)} returns a pointer to a filled 2946 \code{Q00\_PSI\_PHI} structure. 2947 2948 \code{psi} is a pointer to the first set of basis functions 2949 \code{phi} is a pointer to the second set of basis functions; 2950 if both are \nil pointers, nothing is done and the return value 2951 is \nil; if one of the pointers is a \nil pointer, the structure 2952 is initialized using the same set of basis functions for the 2953 first and second set, i.e. \code{phi = psi} or \code{psi = phi} 2954 is used. 2955 2956 \code{quad} is a pointer to a quadrature for the approximation 2957 of the integrals; is \code{quad} is \nil, a quadrature which 2958 is exact of degree \code{psi->degree+phi->degree} is used. 2959 2960 All used \code{Q00\_PSI\_PHI} structures are stored in a linked 2961 list and are identified uniquely by the members \code{psi}, 2962 \code{phi}, and \code{quad}, and for such a combination only 2963 one \code{Q00\_PSI\_PHI} structure is created during runtime; 2964 2965 First, \code{get\_q00\_psi\_phi()} 2966 looks for a matching structure in the linked list; if such a 2967 structure is found a pointer to this structure is returned; 2968 the values are not computed a second time. 2969 Otherwise a new structure is generated, linked to the list, and 2970 the values are computed using the quadrature \code{quad}. 2971\end{descr} 2972 2973\subsection{Data structures and functions for updating coefficient vectors}% 2974\label{S:vector_update}% 2975 2976\sloppypar Besides the general routines \code{update\_real\_vec()}, 2977\code{update\_real\_d\_vec()} and \code{update\_real\_vec\_dow()}, 2978this section presents also easy to use routines for calculation of 2979$L^2$ scalar products between a given function and all basis functions 2980of a finite element space, taken either over the interior of the mesh 2981or over boundary segments. 2982 2983\medskip 2984 2985The following structures hold full information for the assembling 2986of element vectors. They are used by the functions 2987\code{update\_real\_vec()} and \code{update\_real\_d\_vec()} described below. 2988\ddx{EL_VEC_INFO@{\code{EL\_VEC\_INFO}}} 2989\idx{assemblage tools!EL_VEC_INFO@{\code{EL\_VEC\_INFO}}} 2990\ddx{EL_VEC_D_INFO@{\code{EL\_VEC\_D\_INFO}}} 2991\idx{assemblage tools!EL_VEC_D_INFO@{\code{EL\_VEC\_D\_INFO}}} 2992\ddx{EL_VEC_INFO_D@{\code{EL\_VEC\_INFO\_D}}} 2993\idx{assemblage tools!EL_VEC_INFO_D@{\code{EL\_VEC\_INFO\_D}}} 2994\bv\begin{lstlisting}[name=EL_VEC_INFO,label=T:EL_VEC_INFO] 2995typedef struct el_vec_info EL_VEC_INFO; 2996typedef struct el_vec_d_info EL_VEC_D_INFO; 2997typedef struct el_vec_info_d EL_VEC_INFO_D; 2998 2999typedef const EL_REAL_VEC * 3000(*EL_VEC_FCT)(const EL_INFO *el_info, void *fill_info); 3001 3002typedef struct el_vec_info EL_VEC_INFO; 3003struct el_vec_info 3004{ 3005 const FE_SPACE *fe_space; 3006 3007 BNDRY_FLAGS dirichlet_bndry; 3008 REAL factor; 3009 3010 EL_VEC_FCT el_vec_fct; 3011 void *fill_info; 3012 3013 FLAGS fill_flag; 3014}; 3015 3016typedef const EL_REAL_D_VEC * 3017(*EL_VEC_D_FCT)(const EL_INFO *el_info, void *fill_info); 3018 3019typedef struct el_vec_d_info EL_VEC_D_INFO; 3020struct el_vec_d_info 3021{ 3022 const FE_SPACE *fe_space; 3023 3024 BNDRY_FLAGS dirichlet_bndry; 3025 REAL factor; 3026 3027 EL_VEC_D_FCT el_vec_fct; 3028 void *fill_info; 3029 3030 FLAGS fill_flag; 3031}; 3032 3033typedef const EL_REAL_VEC_D * 3034(*EL_VEC_FCT_D)(const EL_INFO *el_info, void *fill_info); 3035 3036typedef struct el_vec_info_d EL_VEC_INFO_D; 3037struct el_vec_info_d 3038{ 3039 const FE_SPACE *fe_space; 3040 3041 BNDRY_FLAGS dirichlet_bndry; 3042 REAL factor; 3043 3044 EL_VEC_FCT_D el_vec_fct; 3045 void *fill_info; 3046 3047 FLAGS fill_flag; 3048}; 3049\end{lstlisting}\ev 3050Description: 3051\begin{descr} 3052\kitem{fe\_space} the underlying finite element space 3053%% 3054\kitem{dirichlet\_bndry} a bit mask marking the boundary segments 3055which are subject to Dirichlet boundary conditions, see also 3056\secref{S:boundary}. 3057%% 3058\kitem{factor} is a multiplier for the element contributions; usually 3059\code{factor} is \code{1} or \code{-1}. 3060%% 3061\kitem{el\_vec\_fct} is a pointer to a function for the computation of 3062the element vector; \code{el\_vec\_fct(el\_info, fill\_info)} returns 3063a pointer to an element vector of the respective data type, see e.g. 3064\verb|EL_REAL_VEC| on page \pageref{T:EL_REAL_VEC}. This vector stores 3065the computed values for the element described by \code{el\_info}. 3066\code{fill\_info} is a pointer to data needed by 3067\code{el\_vec\_fct()}; the function has to provide memory for storing 3068the element vector, which can be overwritten on the next call. 3069%% 3070\kitem{fill\_info} pointer to data needed by \code{el\_vec\_fct()}; is 3071the second argument of this function. 3072%% 3073\kitem{fill\_flag} the flag for the mesh traversal for assembling the 3074vector. 3075\end{descr} 3076The following function does the update of vectors by assembling element 3077contributions during mesh traversal; information for computing 3078the element vectors is held in a \code{EL\_VEC[\_D]\_INFO} structure: 3079\fdx{update_real_vec()@{\code{update\_real\_vec()}}} 3080\idx{assemblage tools!update_real_vec()@{\code{update\_real\_vec()}}} 3081\fdx{update_real_d_vec()@{\code{update\_real\_d\_vec()}}} 3082\idx{assemblage tools!update_real_d_vec()@{\code{update\_real\_d\_vec()}}} 3083\fdx{update_real_vec_dow()@{\code{update\_real\_vec\_dow()}}} 3084\idx{assemblage tools!update_real_vec_dow()@{\code{update\_real\_vec\_dow()}}} 3085\bv\begin{lstlisting} 3086void update_real_vec(DOF_REAL_VEC *dv, const EL_VEC_INFO *vec_info); 3087void update_real_d_vec(DOF_REAL_D_VEC *dv, const EL_VEC_D_INFO *vec_info); 3088void update_real_vec_dow(DOF_REAL_VEC_D *dv, const EL_VEC_INFO_D *vec_info) 3089\end{lstlisting}\ev 3090\begin{descr} 3091 \kitem{update\_real[\_d]\_vec[\_dow](dv, info)} updates the vector 3092 \code{dr} by traversing the underlying mesh and assembling the 3093 element contributions into the vector; information about the 3094 computation of element vectors and connection of local and global 3095 DOFs is stored in \code{info}. 3096 3097 The flags for the mesh traversal of the mesh 3098 \code{dv->fe\_space->mesh} are stored at \code{info->fill\_flags} 3099 which specifies the elements to be visited and information that 3100 should be present on the elements for the calculation of the element 3101 vectors and boundary information (if \code{info->get\_bound} is not 3102 \nil). 3103 3104 On the elements, information about the global DOFs is accessed by 3105 \code{info->get\_dof} using \code{info->admin}; the boundary type of 3106 the DOFs is accessed by \code{info->get\_bound} if 3107 \code{info->get\_bound} is not a \nil pointer; then the element 3108 vector is computed by \code{info->el\_vec\_fct(el\_info, 3109 info->fill\_info)}; these contributions are finally added to 3110 \code{dv} multiplied by \code{info->factor} by a call of 3111 \code{add\_element[\_d]\_vec[\_dow]()} with all information about global 3112 DOFs, the element vector, and boundary types, if available; 3113 3114 \code{update\_real[\_d]\_vec[\_dow]()} only adds element 3115 contributions; this makes several calls for the assemblage of one 3116 vector possible; before the first call, the vector should be set to 3117 zero by a call of \code{dof\_set[\_d|\_dow](0.0, dv)}. 3118\end{descr} 3119 3120\paragraph{$L^2$- and $H^1$-scalar- products over the bulk phase} 3121In many applications, the load vector is just the $L^2$- or 3122$H^1$-scalar-product of a given function with all basis functions of 3123the finite element space or this scalar product is a part of the right 3124hand side. Such a scalar product can be directly assembled by the 3125following functions. 3126 3127\paragraph{Prototypes} 3128\fdx{L2scp_fct_bas()@{\code{L2scp\_fct\_bas()}}} 3129\idx{assemblage tools!L2scp_fct_bas()@{\code{L2scp\_fct\_bas()}}} 3130\fdx{L2scp_fct_bas_d()@{\code{L2scp\_fct\_bas\_d()}}} 3131\idx{assemblage tools!L2scp_fct_bas_d()@{\code{L2scp\_fct\_bas\_d()}}} 3132\fdx{L2scp_fct_bas_dow()@{\code{L2scp\_fct\_bas\_dow()}}} 3133\idx{assemblage tools!L2scp_fct_bas_dow()@{\code{L2scp\_fct\_bas\_dow()}}} 3134\fdx{L2scp_fct_bas_loc()@{\code{L2scp\_fct\_bas\_loc()}}} 3135\idx{assemblage tools!L2scp_fct_bas_loc()@{\code{L2scp\_fct\_bas\_loc()}}} 3136\fdx{L2scp_fct_bas_loc_dow()@{\code{L2scp\_fct\_bas\_loc\_dow()}}} 3137\idx{assemblage tools!L2scp_fct_bas_loc_dow()@{\code{L2scp\_fct\_bas\_loc\_dow()}}} 3138\fdx{H1scp_fct_bas()@{\code{H1scp\_fct\_bas()}}} 3139\idx{assemblage tools!H1scp_fct_bas()@{\code{H1scp\_fct\_bas()}}} 3140\fdx{H1scp_fct_bas_dow()@{\code{H1scp\_fct\_bas\_dow()}}} 3141\idx{assemblage tools!H1scp_fct_bas_dow()@{\code{H1scp\_fct\_bas\_dow()}}} 3142%% 3143\bv\begin{lstlisting} 3144void L2scp_fct_bas(FCT_AT_X f, const QUAD *quad, DOF_REAL_VEC *fh); 3145void L2scp_fct_bas_d(FCT_D_AT_X f, const QUAD *, DOF_REAL_D_VEC *fh); 3146void L2scp_fct_bas_dow(FCT_D_AT_X f, const QUAD *quad, DOF_REAL_VEC_D *fh); 3147 3148void L2scp_fct_bas_loc(DOF_REAL_VEC *fh, 3149 LOC_FCT_AT_QP f, void *f_data, FLAGS fill_flag, 3150 const QUAD *quad); 3151void L2scp_fct_bas_loc_d(DOF_REAL_D_VEC *fh, 3152 LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag, 3153 const QUAD *quad); 3154void L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh, 3155 LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag, 3156 const QUAD *quad); 3157 3158void H1scp_fct_bas(GRD_FCT_AT_X grd_f, 3159 const QUAD *quad, DOF_REAL_VEC *fh); 3160void H1scp_fct_bas_d(GRD_FCT_D_AT_X grd_f, 3161 const QUAD *quad, DOF_REAL_D_VEC *fh); 3162void H1scp_fct_bas_dow(GRD_FCT_D_AT_X grd_f, 3163 const QUAD *quad, DOF_REAL_VEC_D *fh); 3164 3165void H1scp_fct_bas_loc(DOF_REAL_VEC *fh, 3166 GRD_LOC_FCT_AT_QP grd_f, void *fd, 3167 FLAGS fill_flag, const QUAD *quad); 3168void H1scp_fct_bas_loc_d(DOF_REAL_VEC_D *fh, 3169 GRD_LOC_FCT_D_AT_QP grd_f, void *fd, 3170 FLAGS fill_flag, const QUAD *quad); 3171void H1scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh, 3172 GRD_LOC_FCT_D_AT_QP grd_f, void *fd, 3173 FLAGS fill_flag, const QUAD *quad); 3174\end{lstlisting}\ev 3175 3176\paragraph{Descriptions} 3177\begin{descr} 3178 \kitem{L2scp\_fct\_bas(f, quad, fh)} 3179 \kitem{L2scp\_fct\_bas\_d(f, quad, fh)} 3180 \kitem{L2scp\_fct\_bas\_dow(f, quad, fh)} 3181 %% 3182 Approximate the $L^2$ scalar products of a given function with all 3183 basis functions by numerical quadrature and add the corresponding 3184 values to a DOF vector 3185 3186 \code{f} is a pointer for the evaluation of the given 3187 function in world coordinates $x$ and returns the 3188 value of that function at $x$; if \code{f} is a \nil 3189 pointer, nothing is done; 3190 3191 \code{fh} is the DOF vector where at the \code{i}--th entry 3192 the approximation of the $L^2$ scalar product of the given 3193 function with the \code{i}--th global basis function of 3194 \code{fh->fe\_space} is added; 3195 3196 \code{quad} is the quadrature for the approximation of 3197 the integral on each leaf element of \code{fh->fe\_space->mesh}; 3198 if \code{quad} is a \nil pointer, a default quadrature 3199 which is exact of degree \code{2*fh->fe\_space->bas\_fcts->degree-2} 3200 is used. 3201 3202 The integrals are approximated by looping over all 3203 leaf elements, computing the approximations to the element 3204 contributions and adding these values to the vector \code{fh} 3205 by \code{add\_element\_vec()}. 3206 3207 The vector \code{fh} is \emph{not} initialized with \code{0.0}; 3208 only the new contributions are added. 3209 %% 3210\kitem{L2scp\_fct\_bas\_d(fd, quad, fhd)} approximates the $L^2$ scalar 3211 products of a given vector valued function with all scalar valued 3212 basis functions by numerical quadrature and adds the 3213 corresponding values to a vector valued DOF vector; 3214 3215 \code{fd} is a pointer for the evaluation of the given function 3216 in world coordinates $x$; \code{fd(x, fx)} returns a pointer to 3217 a vector storing the value at \code{x}; if \code{fx} is not 3218 \nil, the value is stored at \code{fx} otherwise the function 3219 has to provide memory for storing this vector, which can be 3220 overwritten on the next call; if \code{fd} is a \nil pointer, 3221 nothing is done; 3222 3223 \code{fhd} is the DOF vector where at the \code{i}--th entry (a 3224 \code{REAL\_D} vector) the approximation of the $L^2$ scalar 3225 product of the given vector valued function with the \code{i}--th 3226 global (scalar valued) basis function of \code{fhd->fe\_space} 3227 is added; 3228 3229 \code{quad} is the quadrature for the approximation of 3230 the integral on each leaf element of \code{fhd->fe\_space->mesh}; 3231 if \code{quad} is a \nil pointer, a default quadrature 3232 which is exact of degree \code{2*fhd->fe\_space->bas\_fcts->degree-2} 3233 is used. 3234 3235 The integrals are approximated by looping over all 3236 leaf elements, computing the approximations to the element 3237 contributions and adding these values to the vector \code{fhd} 3238 by \code{add\_element\_d\_vec()}. 3239 3240 The vector \code{fhd} is \emph{not} initialized with 3241 $(\code{0.0},\ldots,\code{0.0})$; only the new contributions are 3242 added. 3243\kitem{L2scp\_fct\_bas\_dow(fd, quad, fhd)} 3244\kitem{L2scp\_fct\_bas\_loc(fh, f\_at\_qp, fct\_data, fill\_flag, quad)} 3245\kitem{L2scp\_fct\_bas\_loc\_dow(fh, f\_at\_qp, ud, fill\_flag, quad)} 3246\kitem{H1scp\_fct\_bas(grd\_f, quad, fh)} 3247\kitem{H1scp\_fct\_bas\_dow(grd\_fd, quad, fhd)} 3248\end{descr} 3249 3250\subsection{Boundary conditions}\label{S:boundary_conditions} 3251%% 3252The following six functions act as a front-end to the functions 3253explained further below, therefore we refer the reader to 3254\secref{S:dirichlet_bound}, \ref{S:neumann_bound} and 3255\ref{S:robin_bound} for a deeper discussion of the implementation of 3256Dirichlet, Neumann and Robin boundary conditions within \ALBERTA. 3257 3258% aus l2scp.c 3259\fdx{boundary_conditions()@{\code{boundary\_conditions()}}} 3260\idx{assemblage tools!boundary_conditions()@{\code{boundary\_conditions()}}} 3261\fdx{boundary_conditions_loc()@{\code{boundary\_conditions\_loc()}}} 3262\idx{assemblage tools!boundary_conditions_loc()@{\code{boundary\_conditions\_loc()}}} 3263\fdx{boundary_conditions_dow()@{\code{boundary\_conditions\_dow()}}} 3264\idx{assemblage tools!boundary_conditions_dow()@{\code{boundary\_conditions\_dow()}}} 3265\fdx{boundary_conditions_loc_dow()@{\code{boundary\_conditions\_loc\_dow()}}} 3266\idx{assemblage tools!boundary_conditions_loc_dow()@{\code{boundary\_conditions\_loc\_dow()}}} 3267\bv\begin{lstlisting} 3268bool boundary_conditions(DOF_MATRIX *matrix, DOF_REAL_VEC *fh, 3269 DOF_REAL_VEC *uh, DOF_SCHAR_VEC *bound, 3270 const BNDRY_FLAGS dirichlet_segment, 3271 REAL (*g)(const REAL_D x), 3272 REAL (*gn)(const REAL_D x, const REAL_D normal), 3273 REAL alpha_r, const WALL_QUAD *wall_quad); 3274bool boundary_conditions_loc(DOF_MATRIX *matrix, DOF_REAL_VEC *fh, 3275 DOF_REAL_VEC *uh, DOF_SCHAR_VEC *bound, 3276 const BNDRY_FLAGS dirichlet_segment, 3277 LOC_FCT_AT_QP g_at_qp, LOC_FCT_AT_QP gn_at_qp, 3278 void *app_data, FLAGS fill_flags, 3279 REAL alpha_r, const WALL_QUAD *wall_quad); 3280bool boundary_conditions_d(DOF_MATRIX *matrix, DOF_REAL_D_VEC *fh, 3281 DOF_REAL_D_VEC *uh, DOF_SCHAR_VEC *bound, 3282 const BNDRY_FLAGS dirichlet_segment, 3283 const REAL *(*g)(const REAL_D x, REAL_D res), 3284 const REAL *(*gn)(const REAL_D x, 3285 const REAL_D normal, 3286 REAL_D res), 3287 REAL alpha_r, const WALL_QUAD *wall_quad); 3288bool boundary_conditions_loc_d(DOF_MATRIX *matrix, DOF_REAL_D_VEC *fh, 3289 DOF_REAL_D_VEC *uh, DOF_SCHAR_VEC *bound, 3290 const BNDRY_FLAGS dirichlet_segment, 3291 LOC_FCT_D_AT_QP g_at_qp, 3292 LOC_FCT_D_AT_QP gn_at_qp, 3293 void *app_data, FLAGS fill_flags, 3294 REAL alpha_r, const WALL_QUAD *wall_quad); 3295bool boundary_conditions_dow(DOF_MATRIX *matrix, DOF_REAL_VEC_D *fh, 3296 DOF_REAL_VEC_D *uh, DOF_SCHAR_VEC *bound, 3297 const BNDRY_FLAGS dirichlet_segment, 3298 const REAL *(*g)(const REAL_D x, REAL_D res), 3299 const REAL *(*gn)(const REAL_D x, 3300 const REAL_D normal, 3301 REAL_D res), 3302 REAL alpha_r, const WALL_QUAD *wall_quad); 3303bool boundary_conditions_loc_dow(DOF_MATRIX *matrix, DOF_REAL_VEC_D *fh, 3304 DOF_REAL_VEC_D *uh, DOF_SCHAR_VEC *bound, 3305 const BNDRY_FLAGS dirichlet_segment, 3306 LOC_FCT_D_AT_QP g_at_qp, 3307 LOC_FCT_D_AT_QP gn_at_qp, 3308 void *app_data, FLAGS fill_flags, 3309 REAL alpha_r, const WALL_QUAD *wall_quad); 3310\end{lstlisting}\ev 3311Description: These ``compound'' functions implement Dirichlet, Neumann 3312or Robin boundary conditions, and optionally perform a mean-value 3313correction of the ``right hand side'' in the context of pure Neumann 3314boundary conditions if $\code{alpha\_r}<0$ (in order to satisfy the 3315conditions for the ``right hand side'' which may be violated in the 3316discrete context because of quadrature errors). 3317 3318For the differences between the code{\dots\_loc()} and 3319non-\code{\dots\_loc()} versions the reader is referred to the section 3320dealing with \code{dirichlet\_bound\_loc()} (see 3321\secref{S:dirichlet_bound}). A brief discussion of the calling 3322convention for the various functions pointers passed to the library 3323functions can also be found in \secref{S:user_fcts}. 3324\begin{description} 3325\item[Parameters]~\hfill 3326 \begin{descr} 3327 \kitem{matrix} As explained in \secref{S:robin_bound}, passed on 3328 to \code{robin\_bound()}. 3329 %% 3330 \kitem{fh} As explained in \secref{S:dirichlet_bound}, 3331 \secref{S:neumann_bound} and \secref{S:robin_bound}. Passed on to 3332 \code{dirichlet\_bound()} and \code{bndry\_L2scp\_fct\_bas()}. 3333 %% 3334 \kitem{uh} As explained in \secref{S:dirichlet_bound}. Passed on 3335 to \code{dirichlet\_bound()}. 3336 %% 3337 \kitem{bound} As explained in \secref{S:dirichlet_bound}. Passed on 3338 to \code{dirichlet\_bound()}. 3339 %% 3340 \kitem{dirichlet\_segment} As explained in 3341 \secref{S:dirichlet_bound}. Passed on to 3342 \code{dirichlet\_bound()}. The respective bit-masks passed to 3343 \code{bndry\_L2scp\_fct\_bas()} and \code{robin\_bound()} are 3344 computed as bit-wise complement of \code{dirichlet\_segment}. See 3345 also \secref{S:boundary}. 3346 %% 3347 \kitem{g} As explained in \secref{S:dirichlet_bound}. Passed on 3348 to \code{dirichlet\_bound()}. 3349 %% 3350 \kitem{gn} As explained in \secref{S:neumann_bound}, 3351 \secref{S:robin_bound}. Passed on to 3352 \code{bndry\_L2scp\_fct\_bas()}. 3353 %% 3354 \kitem{app\_data} \code{\dots\_loc()}-variants only. As explained 3355 in \secref{S:dirichlet_bound}, \secref{S:user_fcts}. Passed on as 3356 application-data pointer to the application provided function 3357 hooks. 3358 %% 3359 \kitem{fill\_flags} \code{\dots\_loc()}-variants only. Additional 3360 fill-flags needed by \code{g()} or \code{gn}. 3361 %% 3362 \kitem{alpha\_r} As explained in \secref{S:robin_bound}. Passed on 3363 to \code{robin\_bound()}. \code{alpha\_r} is also abused to 3364 request an automatic mean-value correction of the load-vector: if 3365 \code{alpha\_r} is negative, and Neumann boundary conditions were 3366 imposed on all boundary segments, then 3367 \code{boundary\_conditions()} will attempt such a mean-value 3368 correction in order to keep fulfill the compatibility condition 3369 for the load-vector in the discrete setting. Of course, if the 3370 differential operator has lower order parts, then the load-vector 3371 need not have mean-value $0$. 3372 3373 Robin boundary conditions will only be assembled if 3374 \code{alpha\_r} is strictly larger than $0$. 3375 %% 3376 \kitem{wall\_quad} As explained in \secref{S:robin_bound} and 3377 \secref{S:neumann_bound}. Passed on to \code{robin\_bound()} and 3378 \code{bndry\_L2scp\_fct\_bas()}. 3379 \end{descr} 3380\item[Return Value]~\hfill 3381 3382 \code{true} if any part of the boundary was subject to Dirichlet 3383 boundary conditions. 3384\end{description} 3385 3386\subsubsection{Dirichlet boundary conditions} 3387\label{S:dirichlet_bound} 3388 3389For the solution of the discrete system \mathref{book:E:LuF} on page 3390\pageref{book:E:LuF} derived in \secref{book:S:Dis2ndorder}, we have 3391to set the Dirichlet boundary values for all Dirichlet DOFs. Usually, 3392we take for the approximation $g_h$ of $g$ the interpolant of $g$, 3393i.e. $g_h = I_h g$ and we have to copy the coefficients of $g_h$ at 3394the Dirichlet DOFs to the start value for an iterative solver. This 3395way the first matrix-vector operation performed by an iterative solver 3396will have the effect to transform an inhomogeneous Dirichlet boundary 3397problem to a homogeneous one by applying the differential operator to 3398the boundary values and subtracting the result from the ``right hand 3399side''. Whether or not it is also necessary to copy these 3400coefficients to the load vector depends on how the matrices act on the 3401coefficients: 3402 3403\begin{itemize} 3404\item If the matrix-rows corresponding to Dirichlet-nodes 3405 $k_1,\,\dots,\,k_M$ have been replaced by unit-vectors $e_{k_l}$ 3406 $(1\leq l\leq M)$, then it is also necessary to copy the Dirichlet 3407 nodes to the load vector (compare \mathref{book:E:Fright} on page 3408 \pageref{book:E:Fright}). Copying the coefficients of $g_h$ at the 3409 Dirichlet DOFs to the initial guess will result in an initial 3410 residual (and then for all subsequent residuals) which is zero at 3411 all Dirichlet \code{DOF}s. 3412 3413 This is the case when Dirichlet bit-masks have been copied to 3414 \code{OPERATOR\_INFO.dirichlet\_bndry} (compare \secref{S:boundary} 3415 and \ref{T:OPERATOR_INFO}); the resulting \code{DOF\_MATRIX} will 3416 then be assembled (\secref{S:matrix_assemblage}) in this way, 3417 replacing any row corresponding to a Dirichlet-node by the 3418 corresponding unit-vector. 3419\item If the matrix-rows corresponding to Dirichlet-nodes have not 3420 been replaced by unit-vectors, then it is still possible to solve a 3421 Dirichlet-problem by passing a \code{DOF\_SCHAR\_VEC} to the 3422 matrix-vector routines (compare \secref{S:solver}, describing the 3423 linear solvers available in \ALBERTA). However, in this case the 3424 matrix-vector subroutines will clear all Dirichlet-nodes to zero, 3425 see \secref{S:DOF_BLAS}. Therefore, in this case it is necessary to 3426 clear the coefficients of the ``right hand side'' which correspond 3427 to Dirichlet-nodes. See the \exampleref{E:CLEARING_DIRICHLET_NODES} 3428 for simple examples how to perform this task. 3429\end{itemize} 3430Note that the matrices generated this way -- either by clearing 3431Dirichlet-rows or by masking out Dirichlet rows -- are not symmetric 3432(compare also \mathref{book:E:matrix} on page \pageref{book:E:matrix}) 3433even if the underlying differential operator is symmetric. However, 3434the restriction of the matrix to the space spanned by the 3435non-Dirichlet \code{DOFs} \emph{is} symmetric, so any iterative solver 3436for symmetric matrices will work, provided one either sets the 3437Dirichlet-values also in the load-vector (if the matrix acts as 3438identity on the Dirichlet \code{DOF}s) or clears the Dirichlet-nodes 3439in the load-vector (if the matrix acts as zero-operator on the 3440Dirichlet \code{DOF}s). 3441 3442The following functions will set Dirichlet boundary values 3443for all DOFs on the Dirichlet boundary, using an interpolation 3444of the boundary values $g$: 3445\fdx{dirichlet_bound()@{\code{dirichlet\_bound()}}} 3446\idx{assemblage tools!dirichlet_bound()@{\code{dirichlet\_bound()}}} 3447\fdx{dirichlet_bound_loc()@{\code{dirichlet\_bound\_loc()}}} 3448\idx{assemblage tools!dirichlet_bound_loc()@{\code{dirichlet\_bound\_loc()}}} 3449\fdx{dirichlet_bound_d()@{\code{dirichlet\_bound\_d()}}} 3450\idx{assemblage tools!dirichlet_bound_d()@{\code{dirichlet\_bound\_d()}}} 3451\fdx{dirichlet_bound_loc_d()@{\code{dirichlet\_bound\_loc\_d()}}} 3452\idx{assemblage tools!dirichlet_bound_loc_d()@{\code{dirichlet\_bound\_loc\_d()}}} 3453\fdx{dirichlet_bound_loc_dow()@{\code{dirichlet\_bound\_loc\_dow()}}} 3454\idx{assemblage tools!dirichlet_bound_loc_dow()@{\code{dirichlet\_bound\_loc\_dow()}}} 3455\bv\begin{lstlisting} 3456bool dirichlet_bound(DOF_REAL_VEC *fh, DOF_REAL_VEC *uh, 3457 DOF_SCHAR_VEC *bound, 3458 const BNDRY_FLAGS dirichlet_segments, 3459 REAL (*g)(const REAL_D)); 3460bool dirichlet_bound_d(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh, 3461 DOF_SCHAR_VEC *bound, 3462 const BNDRY_FLAGS dirichlet_segments, 3463 const REAL *(*g)(const REAL_D, REAL_D)); 3464bool dirichlet_bound_dow(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh, 3465 DOF_SCHAR_VEC *bound, 3466 const BNDRY_FLAGS dirichlet_segments, 3467 const REAL *(*g)(const REAL_D, REAL_D)); 3468bool dirichlet_bound_loc(DOF_REAL_VEC *fh, DOF_REAL_VEC *uh, 3469 DOF_SCHAR_VEC *bound, 3470 const BNDRY_FLAGS dirichlet_segments, 3471 LOC_FCT_AT_QP g, void *ud, FLAGS fill_flags); 3472bool dirichlet_bound_loc_d(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh, 3473 DOF_SCHAR_VEC *bound, 3474 const BNDRY_FLAGS dirichlet_segments, 3475 LOC_FCT_D_AT_QP g, void *ud, 3476 FLAGS fill_flags); 3477bool dirichlet_bound_loc_dow(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh, 3478 DOF_SCHAR_VEC *bound, 3479 const BNDRY_FLAGS dirichlet_segments, 3480 LOC_FCT_D_AT_QP g, void *ud, 3481 FLAGS fill_flags); 3482\end{lstlisting}\ev 3483 3484\paragraph{Descriptions} 3485\begin{descr} 3486 \kitem{dirichlet\_bound(fh, uh, bound, dirichlet\_segments, g)} sets 3487 Dirichlet boundary values for all DOFs on all leaf elements of 3488 \code{fh->fe\_space->mesh} or \code{uh->fe\_space->mesh}; values at 3489 DOFs not belonging to the Dirichlet boundary are not changed by this 3490 function. 3491 3492 Boundary values are set element-wise on the leaf elements. The 3493 boundary type of the walls of an element is accessed through the 3494 function \code{wall\_bound(el\_info, wall)}. If the corresponding 3495 bit is set in \code{dirichlet\_segments}, then the local 3496 interpolation operator of the basis functions underlying 3497 \code{fh/uh->fe\_space} is invoked to compute the coefficients of 3498 the \code{DOF}s located on that wall. 3499 3500 This variant of the \code{dirichlet\_bound...()} is rather 3501 simplistic; the \code{dirichlet\_bound\_loc..()} pass more 3502 information to the function implementing the boundary values and 3503 also allow for manipulating the amount of information available 3504 while looping over the mesh. 3505 %% 3506 \begin{description} 3507 \item[Parameters]~\hfill 3508 \begin{description} 3509 \item[\code{fh}, \code{uh}] are vectors where Dirichlet boundary 3510 values should be set (usually, \code{fh} stores the load vector 3511 and \code{uh} an initial guess for an iterative solver); one of 3512 \code{fh} and \code{uh} may be a \nil pointer; if both arguments 3513 are \nil pointers, nothing is done, except of filling the 3514 \code{bound} argument, it that is non \nil; if both arguments are 3515 not \nil, \code{fh->fe\_space} must equal \code{uh->fe\_space}. 3516 \item[\code{bound}] is a vector for storing the boundary type for 3517 each used DOF; \code{bound} may be \nil; if it is not \nil, the 3518 \code{i}-th entry of the vector is filled with the boundary 3519 type of the \code{i}-th DOF. \code{bound->fe\_space} must be 3520 the same as \code{fh}'s or \code{uh}'s \code{fe\_space}. 3521 \item[\code{dirichlet\_segments}] Bit-mask marking those parts of 3522 the boundary which are subject to Dirichlet boundary 3523 conditions. If bit number $k>0$ is set in 3524 \code{dirichlet\_segments} then the part of the boundary with 3525 boundary classification $k$ is marked as Dirichlet boundary. 3526 Compare \secref{S:boundary}. 3527 \item[\code{REAL (*g)(const REAL\_D arg)}] is a pointer to a 3528 function for the evaluation of the boundary data; if \code{g} 3529 is a \nil pointer, all coefficients at Dirichlet DOFs are set 3530 to \code{0.0}. \code{arg} are the Cartesian co-ordinates where 3531 the value of \code{g} should be computed. 3532 \end{description} 3533 \item[Return Value] \code{true} if any part of the boundary of the 3534 mesh is subject to Dirichlet boundary conditions, as indicated by 3535 the argument \code{dirichlet\_segments}, \code{false} otherwise. 3536 \end{description} 3537 %% 3538 \hrulefill 3539 %% 3540 \kitem{dirichlet\_bound\_d(fh, uh, bound, dirichlet\_segments, g)} 3541 does the same as \code{dirichlet\_bound()}, but \code{fh} and 3542 \code{uh} are \code{DOF\_REAL\_D\_VEC} vectors. 3543 3544 The calling convention for \code{const REAL (*g)(const REAL\_D arg, 3545 REAL\_D result)} is such that \code{g} must allow for 3546 \code{result} being a \nil-pointer. If so, a pointer to a statically 3547 allocated storage area must be returned, otherwise \code{result} 3548 must be filled with the value of \code{g} at the evaluation point 3549 \code{arg}, see \exampleref{E:user_fct_d} in \secref{S:user_fcts}. 3550 %% 3551 Otherwise everything works as for \code{dirichlet\_bound()}, see 3552 above for the documentation. 3553 %% 3554 3555 \hrulefill 3556 %% 3557 \kitem{dirichlet\_bound\_dow(fh, uh, bound, dirichlet\_segments, g)} 3558 does the same as \code{dirichlet\_bound\_d()}, but \code{fh} and 3559 \code{uh} are \code{DOF\_REAL\_VEC\_D} vectors, that is, \code{uh} 3560 and \code{fh} may belong to a finite element space which is a direct 3561 sum, composed of several finite element spaces (note the location of 3562 the \code{\_D} suffix in the data-type names 3563 \code{DOF\_REAL\_VEC\_D} and \code{DOF\_REAL\_D\_VEC}!). The calling 3564 convention for \code{const REAL (*g)(const REAL\_D arg, REAL\_D 3565 result)} is the same as explained above for 3566 \code{dirichlet\_bound\_d()}. 3567 3568 \hrulefill 3569 %% 3570 \kitem{dirichlet\_bound\_loc(fh, uh, bound, dirichlet\_segments, g, 3571 ud, fill\_flags)} 3572 %% 3573 This function differs from its counterpart without the 3574 \code{\_loc}-suffix only in the calling convention for the function 3575 implementing the Dirichlet boundary conditions. We document only the 3576 differing or additional arguments here and refer the reader to the 3577 documentation of \code{dirichlet\_bound()} above: 3578 \begin{description} 3579 \item[\code{REAL (*g)(const EL\_INFO *el\_info, const QUAD *quad, 3580 int iq, void *ud)}] The function pointer to the function 3581 implementing the Dirichlet boundary values. In contrast to the 3582 corresponding function-pointer passed to \code{dirichlet\_bound()} 3583 this function is invoked with a co-dimension $1$ quadrature rule 3584 (compare the \code{interpol}-hooks in the \code{BAS\_FCTS} 3585 structure, \ref{T:BAS_FCTS}, and the definition of the 3586 \code{QUAD}-structure, \ref{T:QUAD}) and a quadrature point, and 3587 first and not least with a filled \code{EL\_INFO}-structure. 3588 3589 This means that \code{g} has full-access to all the information 3590 available through the \code{EL\_INFO} element descriptor. The 3591 amount of data filled-in during mesh-traversal can additionally be 3592 controlled by setting specific fill-flags through the argument 3593 \code{fill\_flags}, which is passed as last argument to 3594 \code{dirichlet\_bound\_loc()}. The last argument \code{ud} to 3595 \code{g} is the same as the pointer \code{ud} passed as pre-last 3596 argument to \code{dirichlet\_bound\_loc()} and may be used by an 3597 application to reduce the amount of global variables and thus the 3598 potential of bugs implied by the use of global variables. The 3599 following simple example shows how to get hold of the Cartesian 3600 co-ordinates of the quadrature point, and how to use, e.g., the 3601 boundary classification available through the \code{EL\_INFO} 3602 structure: 3603 %% 3604 \begin{example} 3605 \label{E:loc_fct_at_qp} 3606 \bv\begin{lstlisting}[name=DIRICHLET_BOUND_LOC_EXAMPLE,label=C:DIRICHLET_BOUND_LOC_EXAMPLRE] 3607 struct g_data 3608 { 3609 BNDRY_TYPE special_wall_type; /* other stuff */ 3610 }; 3611 3612 REAL g_impl(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud) 3613 { 3614 struct g_data *data = (struct g_data *)ud; 3615 BNDRY_TYPE btype = wall_bound(el_info, quad->sub_splx); 3616 REAL result; 3617 const QUAD_EL_CACHE *qelc = 3618 fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD); 3619 3620 if (btype == data->special_wall_type) { 3621 ... /* do special things */ 3622 return sin(qelc->world[iq][0]; 3623 } else { 3624 ... /* do "normal" things */ 3625 return sin(qelc->world[iq][1]; 3626 } 3627 } 3628 3629 ... /* 1.000.000 lines of code later ... */ 3630 struct g_data g_data_instance = { 42 }; 3631 dirichlet_bound_loc(fh, uh, bound, dirichlet_bits, g_impl, &g_data_instance, FILL_COORDS|FILL_MACRO_WALLS); 3632 \end{lstlisting}\ev 3633 \end{example} 3634 %% 3635 \item[\code{ud}] Application-data-pointer, forwarded as \code{ud} 3636 argument to the application supplied \code{g} function-pointer. 3637 %% 3638 \item[\code{fill\_flags}] Additional fill-flags for the loop over 3639 the mesh. The application must make sure that \code{fill\_flags} 3640 contains all flags corresponding to information needed by the 3641 function \code{g()}. 3642 \end{description} 3643 %% 3644 \hrulefill 3645 %% 3646 \kitem{dirichlet\_bound\_loc\_dow(} 3647 \kitem{~~~~~~~~~~~~fh, uh, bound, dirichlet\_segments, g, ud, fill\_flags)} 3648 %% 3649 \kitem{dirichlet\_bound\_loc(fh, uh, bound, dirichlet\_segments, g, ud, fill\_flags)} 3650 %% 3651 These two function differ from \code{dirichlet\_bound\_loc()} only 3652 in the calling convention for 3653 \bv\begin{lstlisting} 3654 const REAL *(*g)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int iq, void *ud). 3655 \end{lstlisting}\ev 3656 As in the example \ref{C:FCT_D_AT_X_ABI} the implementation of 3657 \code{g()} must allow for \code{result} being \nil, returning a 3658 pointer to a static storage area in this case. 3659\end{descr} 3660 3661\begin{example} 3662 \label{E:CLEARING_DIRICHLET_NODES} 3663 This example demonstrates how to clear the Dirichlet-nodes in the 3664 load-vector if Dirichlet boundary conditions are implemented using a 3665 \code{DOF\_SCHAR\_VEC} to mask-out Dirichlet nodes. Note that this 3666 example applies \emph{only} if the \code{DOF\_SCHAR\_VEC} is also 3667 passed to the linear solvers. Otherwise Dirichlet boundary 3668 conditions have to be incorporated into the matrix 3669 \begin{description} 3670 \item[\emph{scalar problem}] 3671 \bv\begin{lstlisting}[name=CLEARING_SCALAR_DIRICHLET_NODES,label=C:CLEARING_SCALAR_DIRICHLET_NODES] 3672 extern REAL g(const REAL_D x); /* defined somewhere else, e.g. */ 3673 extern DOF_REAL_VEC *uh, *fh; /* defined somewhere else, e.g. */ 3674 DOF_SCHAR_VEC *bound = 3675 get_dof_schar_vec("dirichlet mask vector", fh->fe_space); 3676 BNDRY_FLAGS dirichlet_bits; 3677 BNDRY_FLAGS_INIT(dirichlet_bits); 3678 BNDRY_FLAGS_SET(dirichlet_bits, 3); /* e.g. */ 3679 3680 ... /* other stuff */ 3681 3682 dirichlet_bound(NULL, uh, bound, dirichlet_bits, g); 3683 FOR_ALL_DOFS(fh->fe_space->admin, 3684 if (bound->vec[dof] >= DIRICHLET) { 3685 fh->vec[dof] = 0.0; 3686 }); 3687 3688 ... /* other stuff */ 3689 3690 oem_solve_s(matrix, bound, fh, uh, ... /* other parameters */); 3691 \end{lstlisting}\ev 3692 \item[\emph{simple vector valued problem}] 3693 \bv\begin{lstlisting}[name=CLEARING_VECTOR_DIRICHLET_NODES,label=C:CLEARING_VECTOR_DIRICHLET_NODES] 3694 extern const REAL *g(const REAL_D x, REAL_D result); /* defined somewhere else, e.g. */ 3695 extern DOF_REAL_D_VEC *uh, *fh; /* defined somewhere else, e.g. */ 3696 extern DOF_SCHAR_VEC *bound; 3697 extern BNDRY_FLAGS dirichlet_bits; 3698 3699 ... /* other stuff */ 3700 3701 dirichlet_bound_d(NULL, uh, bound, dirichlet_bits, g); 3702 FOR_ALL_DOFS(fh->fe_space->admin, 3703 if (bound->vec[dof] >= DIRICHLET) { 3704 SET_DOW(0.0, fh->vec[dof]); 3705 }); 3706 3707 ... /* other stuff */ 3708 3709 oem_solve_d(matrix, bound, fh, uh, ... /* other parameters */); 3710 \end{lstlisting}\ev 3711 \item[\emph{vector valued problem, using an FE-space which is a direct sum}] 3712 ~\hfill 3713 3714 %% 3715 Note the difference between a \code{DOF\_REAL\_D\_VEC} which 3716 contains \code{DIM\_OF\_WORLD}-sized \code{REAL\_D} vectors and a 3717 \code{DOF\_REAL\_VEC\_D} which contains scalars of type 3718 \code{REAL} if the underlying basis function are themselves 3719 vector-valued, or \code{REAL\_D}-vectors if the underlying basis 3720 functions are scalar-valued. The first code-block of the 3721 \code{FOREACH\_DOF\_DOW}-macro is for the case where the basis 3722 functions are vector-valued (and hence the coefficients are 3723 scalars) and the second code-block is for the case where the basis 3724 functions are scalar-valued (and hence the coefficients are 3725 vectors). The \code{FOREACH\_DOF\_DOW()} macro is further 3726 explained in \secref{S:chain_loops}. 3727 %% 3728 \bv\begin{lstlisting}[name=CLEARING_CHAINED_DIRICHLET_NODES,label=C:CLEARING_CHAINED_DIRICHLET_NODES] 3729 extern const REAL *g(const REAL_D x, REAL_D result); /* defined somewhere else, e.g. */ 3730 extern DOF_REAL_VEC_D *uh, *fh; /* defined somewhere else, e.g. */ 3731 extern DOF_SCHAR_VEC *bound; 3732 extern BNDRY_FLAGS dirichlet_bits; 3733 3734 ... /* other stuff */ 3735 3736 dirichlet_bound_dow(NULL, uh, bound, dirichlet_bits, g); 3737 FOREACH_DOF_DOW(fh->fe_space, 3738 if (bound->vec[dof] >= DIRICHLET) { 3739 fh->vec[dof] = 0.0; 3740 }, 3741 if (bound->vec[dof] >= DIRICHLET) { 3742 SET_DOW(0.0, ((DOF_REAL_D_VEC *)fh)->vec[dof]); 3743 }, 3744 fh = CHAIN_NEXT(fh, DOF_REAL_VEC_D); 3745 bound = CHAIN_NEXT(bound, DOF_SCHAR_VEC)); 3746 3747 ... /* other stuff */ 3748 3749 oem_solve_dow(matrix, bound, fh, uh, ... /* other parameters */); 3750 \end{lstlisting}\ev 3751 \end{description} 3752\end{example} 3753 3754\subsubsection{Neumann boundary conditions} 3755\label{S:neumann_bound} 3756 3757For the implementation of inhomogeneous Neumann boundary conditions it 3758is necessary to compute $L^2$ scalar products between the 3759inhomogeneity and the basis functions on the Neumann boundary 3760segments. The following functions compute the $L^2$ scalar product 3761over the boundary of the domain. They return \code{true} if not all 3762boundary segments of the mesh belong to the segment specified by 3763\code{bndry\_seg}. If \code{bndry\_seg == NULL} then the scalar 3764product is computed over the entire boundary (i.e. over all walls 3765without neighbour). 3766%% 3767Besides computing the $L^2$-scalar product over boundary segments 3768there are also functions to compute the $L^2$-scalar-product over trace 3769meshes (or ``sub-meshes'', see \secref{S:submesh_implementation}). 3770%% 3771For the calling conventions for the application provided function 3772pointers the reader is referred to \secref{S:user_fcts}, and the 3773relevant part of the discussion of \code{dirichlet\_bound\_loc()} in 3774\secref{S:dirichlet_bound}. 3775 3776All function work additive, the contributions of the per-element 3777integrals are added to any data already present in \code{fh}. 3778%% 3779\paragraph{Prototypes} 3780\fdx{bndry_L2scp_fct_bas()@{\code{bndry\_L2scp\_fct\_bas()}}} 3781\idx{assemblage tools!bndry_L2scp_fct_bas()@{\code{bndry\_L2scp\_fct\_bas()}}} 3782\fdx{bndry_L2scp_fct_bas_loc()@{\code{bndry\_L2scp\_fct\_bas\_loc()}}} 3783\idx{assemblage tools!bndry_L2scp_fct_bas_loc()@{\code{bndry\_L2scp\_fct\_bas\_loc()}}} 3784\fdx{bndry_L2scp_fct_bas_dow()@{\code{bndry\_L2scp\_fct\_bas\_dow()}}} 3785\idx{assemblage tools!bndry_L2scp_fct_bas_dow()@{\code{bndry\_L2scp\_fct\_bas\_dow()}}} 3786\fdx{bndry_L2scp_fct_bas_loc_dow()@{\code{bndry\_L2scp\_fct\_bas\_loc\_dow()}}} 3787\idx{assemblage tools!bndry_L2scp_fct_bas_loc_dow()@{\code{bndry\_L2scp\_fct\_bas\_loc\_dow()}}} 3788\fdx{trace_L2scp_fct_bas()@{\code{trace\_L2scp\_fct\_bas()}}} 3789\idx{evaluation tools!trace_L2scp_fct_bas()@{\code{trace\_L2scp\_fct\_bas()}}} 3790\fdx{trace_L2scp_fct_bas_loc()@{\code{trace\_L2scp\_fct\_bas\_loc()}}} 3791\idx{evaluation tools!trace_L2scp_fct_bas_loc()@{\code{trace\_L2scp\_fct\_bas\_loc()}}} 3792\fdx{trace_L2scp_fct_bas_dow()@{\code{trace\_L2scp\_fct\_bas\_dow()}}} 3793\idx{evaluation tools!trace_L2scp_fct_bas_dow()@{\code{trace\_L2scp\_fct\_bas\_dow()}}} 3794\fdx{trace_L2scp_fct_bas_loc_dow()@{\code{trace\_L2scp\_fct\_bas\_loc\_dow()}}} 3795\idx{evaluation tools!trace_L2scp_fct_bas_loc_dow()@{\code{trace\_L2scp\_fct\_bas\_loc\_dow()}}} 3796%% 3797\bv\begin{lstlisting}[label=F:BNDRY_L2SCP_PROTO] 3798bool bndry_L2scp_fct_bas_loc(DOF_REAL_VEC *fh, 3799 LOC_FCT_AT_QP f_at_qp, void *ud, FLAGS fill_flag, 3800 const BNDRY_FLAGS bndry_seg, 3801 const WALL_QUAD *quad); 3802bool bndry_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh, LOC_FCT_D_AT_QP f_at_qp, 3803 void *ud, FLAGS fill_flag, 3804 const BNDRY_FLAGS bndry_seg, 3805 const WALL_QUAD *quad); 3806bool bndry_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh, 3807 const REAL *(*f)(const REAL_D x, 3808 const REAL_D normal, 3809 REAL_D result), 3810 const BNDRY_FLAGS bndry_seg, 3811 const WALL_QUAD *quad); 3812bool bndry_L2scp_fct_bas(DOF_REAL_VEC *fh, 3813 REAL (*f)(const REAL_D x, const REAL_D normal), 3814 const BNDRY_FLAGS bndry_seg, const WALL_QUAD *quad); 3815 3816void trace_L2scp_fct_bas(DOF_REAL_VEC *fh, FCT_AT_X f, 3817 MESH *trace_mesh, const QUAD *quad); 3818void trace_L2scp_fct_bas_loc(DOF_REAL_VEC *fh, 3819 LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag, 3820 MESH *trace_mesh, 3821 const QUAD *quad); 3822void trace_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh, FCT_D_AT_X f, 3823 MESH *trace_mesh, 3824 const QUAD *quad); 3825void trace_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh, 3826 LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag, 3827 MESH *trace_mesh, 3828 const QUAD *quad); 3829\end{lstlisting}\ev 3830 3831\paragraph{Descriptions} 3832 3833\begin{descr} 3834 \kitem{bndry\_L2scp\_fct\_bas()}~\hfill 3835 \begin{description} 3836 \item[Parameters]~\hfill 3837 \kitem{fh} The load-vector to add the boundary integrals to. 3838 %% 3839 \kitem{f} Application supplied ``right hand side''. 3840 %% 3841 \kitem{ud} Data pointer for \code{f} for the 3842 \code{\dots\_loc()}-variants. 3843 %% 3844 \kitem{fill\_flags} Additional fill-flags needed by \code{f} for the 3845 \code{\dots\_loc()}-variants. 3846 %% 3847 \kitem{bndry\_seg} A bit-mask specifying the part of the boundary 3848 which is the domain of integration. See \secref{S:boundary}. 3849 %% 3850 \kitem{quad} Optional application supplied quadrature rule. Maybe 3851 \nil, in which case a default quadrature rule is used. See 3852 \secref{S:WALL_QUAD} for how to obtain such a structure, function 3853 \code{get\_wall\_quad()}. 3854 \item[Return Value]~\hfill 3855 3856 \code{true} if part of the boundary did 3857 \emph{not} belong to the domain of integration. 3858 \end{description} 3859 3860 \hrulefill 3861 3862 \kitem{trace\_L2scp\_fct\_bas()}~\hfill 3863 \begin{description} 3864 \item[Parameters]~\hfill 3865 \begin{descr} 3866 \kitem{fh} The load-vector. 3867 %% 3868 \kitem{f} The user-supplied inhomogeneity. 3869 %% 3870 \kitem{fd} The application data pointer passed on to \code{f} 3871 for the \code{\dots\_loc()}-variants. 3872 %% 3873 \kitem{fill\_flags} Additional fill-flags for the 3874 \code{\dots\_loc()}-variants. 3875 %% 3876 \kitem{trace\_mesh} The domain of integration. 3877 %% 3878 \kitem{quad} A user supplied quadrature rule. May be \nil in 3879 which case a default quadrature rule will be used. The 3880 quadrature rule must have the dimension of \code{trace\_mesh}, 3881 naturally. 3882 %% 3883 \end{descr} 3884 \item[Return Value]~\hfill 3885 \end{description} 3886\end{descr} 3887%Description: 3888%\begin{descr} 3889%\kitem{bndry\_L2scp\_fct\_bas\_loc(fh, f\_at\_qp, ud, fill\_flag, bndry\_seg, quad)} 3890%\fdx{bndry_L2scp_fct_bas_loc()@{\code{bndry\_L2scp\_fct\_bas\_loc()}}} 3891%\idx{assemblage tools!bndry_L2scp_fct_bas_loc()@{\code{bndry\_L2scp\_fct\_bas\_loc()}}} 3892%%% 3893%\kitem{bndry\_L2scp\_fct\_bas\_loc\_dow(fh, f\_at\_qp, ud, fill\_flag, bndry\_seg, quad)} 3894%\fdx{bndry_L2scp_fct_bas_loc_dow()@{\code{bndry\_L2scp\_fct\_bas\_loc\_dow()}}} 3895%\idx{assemblage tools!bndry_L2scp_fct_bas_loc_dow()@{\code{bndry\_L2scp\_fct\_bas\_loc\_dow()}}} 3896 3897 3898%\kitem{bndry\_L2scp\_fct\_bas\_dow(fh, f, bndry\_seg, quad)} 3899%\fdx{bndry_L2scp_fct_bas_dow()@{\code{bndry\_L2scp\_fct\_bas\_dow()}}} 3900%\idx{assemblage tools!bndry_L2scp_fct_bas_dow()@{\code{bndry\_L2scp\_fct\_bas\_dow()}}} 3901%%% 3902%\kitem{bndry\_L2scp\_fct\_bas(fh, f, bndry\_seg, quad)} 3903%\fdx{bndry_L2scp_fct_bas()@{\code{bndry\_L2scp\_fct\_bas()}}} 3904%\idx{assemblage tools!bndry_L2scp_fct_bas()@{\code{bndry\_L2scp\_fct\_bas()}}} 3905 3906%\kitem{trace\_L2scp\_fct\_bas(fh, f, trace\_mesh, quad)} 3907%\fdx{trace_L2scp_fct_bas()@{\code{trace\_L2scp\_fct\_bas()}}} 3908%\idx{evaluation tools!trace_L2scp_fct_bas()@{\code{trace\_L2scp\_fct\_bas()}}} 3909 3910%\kitem{trace\_L2scp\_fct\_bas\_loc(fh, f, fd, fill\_flag, trace\_mesh, quad)} 3911%\fdx{trace_L2scp_fct_bas_loc()@{\code{trace\_L2scp\_fct\_bas\_loc()}}} 3912%\idx{evaluation tools!trace_L2scp_fct_bas_loc()@{\code{trace\_L2scp\_fct\_bas\_loc()}}} 3913 3914%\kitem{trace\_L2scp\_fct\_bas\_dow(fh, f, trace\_mesh, quad)} 3915%\fdx{trace_L2scp_fct_bas_dow()@{\code{trace\_L2scp\_fct\_bas\_dow()}}} 3916%\idx{evaluation tools!trace_L2scp_fct_bas_dow()@{\code{trace\_L2scp\_fct\_bas\_dow()}}} 3917 3918%\kitem{trace\_L2scp\_fct\_bas\_loc\_dow(fh, f, fd, fill\_flag, trace\_mesh, quad)} 3919%\fdx{trace_L2scp_fct_bas_loc_dow()@{\code{trace\_L2scp\_fct\_bas\_loc\_dow()}}} 3920%\idx{evaluation tools!trace_L2scp_fct_bas_loc_dow()@{\code{trace\_L2scp\_fct\_bas\_loc\_dow()}}} 3921%\end{descr} 3922 3923\subsubsection{Robin boundary conditions} 3924\label{S:robin_bound} 3925 3926\fdx{robin_bound()@{\code{robin\_bound()}}} 3927\idx{assemblage tools!robin_bound()@{\code{robin\_bound()}}} 3928\bv\begin{lstlisting} 3929void robin_bound(DOF_MATRIX *matrix, const BNDRY_FLAGS robin_seg, 3930 REAL alpha_r, const WALL_QUAD *wall_quad, 3931 REAL exponent); 3932\end{lstlisting}\ev 3933Description: Incorporate so-called ``Robin boundary'' conditions into 3934the matrix, i.e. a boundary condition of the form 3935\[ 3936\frac{\partial u}{\partial\nu}(x) + \alpha(x)\,u(x) = g(x)\quad\text{ on }\partial\Omega. 3937\] 3938In the context of weak formulations for elliptic second-order PDEs, 3939this results into two boundary integrals, one has to be added to the 3940linear form on the ``right hand side'', and the other one is a 3941contribution to bilinear-form on the ``left hand side'', namely 3942\[ 3943\int_{\partial\Omega} \alpha\,u\,\phi\,do 3944\quad\text{ and }\quad 3945\int_{\partial\Omega}g\,\phi\,do 3946\] 3947The contribution to the right hand side can be assembled using one of 3948the \code{bndry\_L2scp\_fct\_bas()}-variants, the contribution the 3949left hand side should be added to the system matrix. 3950\code{robin\_bound()} implements this for the case where $\alpha$ is 3951actually just a constant coefficient. 3952\begin{descr} 3953\kitem{robin\_bound(matrix, robin\_seg, alpha\_r, wall\_quad, exponent)}~\hfill 3954\begin{description} 3955\item[Parameters]~\hfill 3956 \begin{descr} 3957 \kitem{matrix} The system matrix, the contributions from the Robin 3958 boundary condition are added to \code{matrix}. 3959 %% 3960 \kitem{robin\_segment} A boundary bit-mask, marking all boundary 3961 segments which are subject to the Robin boundary condition. The 3962 position of the bits set in \code{robin\_segment} correspond to 3963 the boundary numbers assigned to the mesh boundary in the macro 3964 triangulation, compare \secref{S:macro_tria} and 3965 \secref{S:boundary}. 3966 %% 3967 \kitem{alpha\_r} The constant coefficient from the Robin boundary 3968 condition. 3969 %% 3970 \kitem{wall\_quad} Optional. A collection of co-dimension $1$ 3971 quadrature formulas for doing the integration. If \code{wall\_quad 3972 == \nil}, then \code{robin\_bound()} chooses a default 3973 quadrature formula, based on the polynomial degree of the 3974 underlying basis-functions. 3975 %% 3976 \kitem{exponent} If \code{exponent > 0.0}, then the boundary 3977 integral will be weighted by the factor $h(T)^{-\code{exponent}}$, 3978 where $h(T)$ denotes the local mesh-width. 3979 \end{descr} 3980\end{description} 3981\end{descr} 3982 3983\subsection{Interpolation into finite element spaces} 3984\label{S:I_FES} 3985 3986In time dependent problems, usually the ``solve'' step in the adaptive 3987method for the adaptation of the initial grid is an interpolation of 3988initial data to the finite element space, i.e. a DOF vector is filled 3989with the coefficient of the interpolant. The following functions are 3990implemented for this task: 3991%% 3992\fdx{interpol()@{\code{interpol()}}} 3993\idx{assemblage tools!interpol()@{\code{interpol()}}} 3994\fdx{interpol_d()@{\code{interpol\_d()}}} 3995\idx{assemblage tools!interpol_d()@{\code{interpol\_d()}}} 3996\fdx{interpol_dow()@{\code{interpol\_dow()}}} 3997\idx{assemblage tools!interpol_dow()@{\code{interpol\_dow()}}} 3998%% 3999\fdx{interpol_loc()@{\code{interpol\_loc()}}} 4000\idx{assemblage tools!interpol()@{\code{interpol\_loc()}}} 4001\fdx{interpol_loc_d()@{\code{interpol\_loc\_d()}}} 4002\idx{assemblage tools!interpol_loc_d()@{\code{interpol\_loc\_d()}}} 4003\fdx{interpol_loc_dow()@{\code{interpol\_loc\_dow()}}} 4004\idx{assemblage tools!interpol_loc_dow()@{\code{interpol\_loc\_dow()}}} 4005%% 4006\bv\begin{lstlisting} 4007void interpol(FCT_AT_X f, DOF_REAL_VEC *fh); 4008void interpol_d(const REAL *(*f)(const REAL_D, REAL_D), DOF_REAL_D_VEC *fh); 4009void interpol_dow(FCT_D_AT_X f, DOF_REAL_VEC_D *fh); 4010void interpol_loc(DOF_REAL_VEC *fh, 4011 LOC_FCT_AT_QP f, void *f_data, FLAGS fill_flags); 4012void interpol_loc_d(DOF_REAL_D_VEC *fh, 4013 LOC_FCT_D_AT_QP f, void *f_data, FLAGS fill_flags); 4014void interpol_loc_dow(DOF_REAL_VEC_D *fh, 4015 LOC_FCT_D_AT_QP f, void *f_data, FLAGS fill_flags); 4016\end{lstlisting}\ev 4017Description: 4018\begin{descr} 4019 \kitem{interpol(f, fh)} computes the coefficients of the interpolant 4020 of a function and stores these in a DOF vector; 4021 4022 Interpolation is done element--wise on the leaf elements of 4023 \code{fh->fe\_space->mesh}; the element interpolation is done by the 4024 function \code{fh->fe\_space->bas\_fcts->interpol()}; the 4025 \code{fill\_flag} of the mesh traversal is 4026 \code{CALL\_LEAF\_EL|FILL\_COORDS} and the function \code{f} must 4027 fit to the needs of \code{fh->fe\_space->bas\_fcts->interpol()}; for 4028 Lagrange elements, \code{(*f)()} is evaluated for all Lagrange nodes 4029 on the element and has to return the value at these nodes (compare 4030 \secref{S:basfct_data}). 4031 4032 \begin{description} 4033 \item[Parameters]~\hfill 4034 \begin{descr} 4035 \kitem{f} is a pointer to a function for the 4036 evaluation of the function to be interpolated; if \code{f} is a \nil 4037 pointer, all coefficients are set to \code{0.0}. 4038 %% 4039 \kitem{fh} is a DOF vector for storing the coefficients; if 4040 \code{fh} is a \nil pointer, nothing is done. 4041 \end{descr} 4042 \end{description} 4043 %% 4044 \kitem{interpol\_d(fd, fhd)} computes the coefficients of the 4045 interpolant of a vector valued function and stores these in a DOF 4046 vector. This version is for the case where the underlying 4047 basis-functions are themselves scalars, consequently the coefficient 4048 vector \code{fh} has the type \code{DOF\_REAL\_D\_VEC}. Otherwise 4049 this function differs from the scalar counter-part only in the 4050 calling convention for the application supplied function \code{f}, 4051 which is the same as for \code{dirichlet\_bound\_d()}, see also 4052 \exampleref{E:user_fct_d} 4053 %% 4054 \kitem{interpol\_dow(fct, uh)} same as \code{interpol\_d()}, but for 4055 the case where the underlying basis function are either scalar- or 4056 \DOW-valued and the finite-element space may optionally be a direct 4057 sum of finite element spaces. 4058 %% 4059 \kitem{interpol\_loc(vec, fct\_at\_qp, app\_data, fill\_flags)} 4060 %% 4061 \kitem{interpol\_loc\_d(vec, fct\_at\_qp, app\_data, fill\_flags)} 4062 %% 4063 \kitem{interpol\_loc\_dow(vec, fct\_at\_qp, app\_data, fill\_flags)} 4064 %% 4065 The \code{\dots\_loc\dots}-variants differ from the other 4066 \code{interpol()}-flavours only in the calling convention for the 4067 application supplied function and the additional \code{fill\_flags} 4068 argument. This has already be explained above, see also 4069 \exampleref{E:loc_fct_at_qp}. 4070\end{descr} 4071\idx{assemblage tools|)} 4072 4073%%% Local Variables: 4074%%% mode: latex 4075%%% TeX-master: "alberta-man" 4076%%% End: 4077