// global right hand side to zero, before we fill them:
cell_matrix = 0;
cell_rhs = 0;
-
- // Then finally assemble the matrix: For the Laplace problem, the matrix
- // on each cell is the integral over the gradients of shape function i
- // and j. Since we do not integrate, but rather use quadrature, this is
- // the sum over all quadrature points of the integrands times the
- // determinant of the Jacobian matrix at the quadrature point times the
- // weight of this quadrature point. You can get the gradient of shape
- // function $i$ at quadrature point q_point by using
- // <code>fe_values.shape_grad(i,q_point)</code>; this gradient is a
- // 2-dimensional vector (in fact it is of type Tensor@<1,dim@>, with
- // here dim=2) and the product of two such vectors is the scalar
- // product, i.e. the product of the two shape_grad function calls is the
- // dot product. This is in turn multiplied by the Jacobian determinant
- // and the quadrature point weight (that one gets together by the call
- // to FEValues::JxW() ). Finally, this is repeated for all shape
- // functions $i$ and $j$:
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
- fe_values.shape_grad (j, q_point) *
- fe_values.JxW (q_point));
-
- // We then do the same thing for the right hand side. Here, the integral
- // is over the shape function i times the right hand side function,
- // which we choose to be the function with constant value one (more
- // interesting examples will be considered in the following
- // programs). Again, we compute the integral by quadrature, which
- // transforms the integral to a sum over all quadrature points of the
- // value of the shape function at that point times the right hand side
- // function, here the constant function equal to one, times the Jacobian
- // determinant times the weight of that quadrature point:
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- cell_rhs(i) += (fe_values.shape_value (i, q_point) *
- 1 *
- fe_values.JxW (q_point));
-
+
+ // Now it is time to start integration over the cell, which we
+ // do by looping over all quadrature points, which we will
+ // number by q_index.
+ for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
+ {
+ // First assemble the matrix: For the Laplace problem, the
+ // matrix on each cell is the integral over the gradients of
+ // shape function i and j. Since we do not integrate, but
+ // rather use quadrature, this is the sum over all
+ // quadrature points of the integrands times the determinant
+ // of the Jacobian matrix at the quadrature point times the
+ // weight of this quadrature point. You can get the gradient
+ // of shape function $i$ at quadrature point with number q_index by
+ // using <code>fe_values.shape_grad(i,q_index)</code>; this
+ // gradient is a 2-dimensional vector (in fact it is of type
+ // Tensor@<1,dim@>, with here dim=2) and the product of two
+ // such vectors is the scalar product, i.e. the product of
+ // the two shape_grad function calls is the dot
+ // product. This is in turn multiplied by the Jacobian
+ // determinant and the quadrature point weight (that one
+ // gets together by the call to FEValues::JxW() ). Finally,
+ // this is repeated for all shape functions $i$ and $j$:
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (fe_values.shape_grad (i, q_index) *
+ fe_values.shape_grad (j, q_index) *
+ fe_values.JxW (q_index));
+
+ // We then do the same thing for the right hand side. Here,
+ // the integral is over the shape function i times the right
+ // hand side function, which we choose to be the function
+ // with constant value one (more interesting examples will
+ // be considered in the following programs).
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_rhs(i) += (fe_values.shape_value (i, q_index) *
+ 1 *
+ fe_values.JxW (q_index));
+ }
// Now that we have the contribution of this cell, we have to transfer
// it to the global matrix and right hand side. To this end, we first
// have to find out which global numbers the degrees of freedom on this
// difference to how we did things in step-3: Instead of using a
// constant right hand side with value 1, we use the object representing
// the right hand side and evaluate it at the quadrature points:
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
- fe_values.shape_grad (j, q_point) *
- fe_values.JxW (q_point));
+ cell_matrix(i,j) += (fe_values.shape_grad (i, q_index) *
+ fe_values.shape_grad (j, q_index) *
+ fe_values.JxW (q_index));
- cell_rhs(i) += (fe_values.shape_value (i, q_point) *
- right_hand_side.value (fe_values.quadrature_point (q_point)) *
- fe_values.JxW (q_point));
+ cell_rhs(i) += (fe_values.shape_value (i, q_index) *
+ right_hand_side.value (fe_values.quadrature_point (q_index)) *
+ fe_values.JxW (q_index));
}
// As a final remark to these loops: when we assemble the local
// contributions into <code>cell_matrix(i,j)</code>, we have to multiply
- // the gradients of shape functions $i$ and $j$ at point q_point and
+ // the gradients of shape functions $i$ and $j$ at point number
+ // q_index and
// multiply it with the scalar weights JxW. This is what actually
- // happens: <code>fe_values.shape_grad(i,q_point)</code> returns a
+ // happens: <code>fe_values.shape_grad(i,q_index)</code> returns a
// <code>dim</code> dimensional vector, represented by a
// <code>Tensor@<1,dim@></code> object, and the operator* that
// multiplies it with the result of
- // <code>fe_values.shape_grad(j,q_point)</code> makes sure that the
+ // <code>fe_values.shape_grad(j,q_index)</code> makes sure that the
// <code>dim</code> components of the two vectors are properly
// contracted, and the result is a scalar floating point number that
// then is multiplied with the weights. Internally, this operator* makes
coefficient.value_list (fe_values.get_quadrature_points(),
coefficient_values);
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int q_index=0; q_index<n_q_indexs; ++q_index)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (coefficient_values[q_point] *
- fe_values.shape_grad(i,q_point) *
- fe_values.shape_grad(j,q_point) *
- fe_values.JxW(q_point));
+ cell_matrix(i,j) += (coefficient_values[q_index] *
+ fe_values.shape_grad(i,q_index) *
+ fe_values.shape_grad(j,q_index) *
+ fe_values.JxW(q_index));
- cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ cell_rhs(i) += (fe_values.shape_value(i,q_index) *
1.0 *
- fe_values.JxW(q_point));
+ fe_values.JxW(q_index));
}
coefficient.value_list (fe_values.get_quadrature_points(),
coefficient_values);
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int q_index=0; q_index<n_q_indexs; ++q_index)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (coefficient_values[q_point] *
- fe_values.shape_grad(i,q_point) *
- fe_values.shape_grad(j,q_point) *
- fe_values.JxW(q_point));
+ cell_matrix(i,j) += (coefficient_values[q_index] *
+ fe_values.shape_grad(i,q_index) *
+ fe_values.shape_grad(j,q_index) *
+ fe_values.JxW(q_index));
- cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ cell_rhs(i) += (fe_values.shape_value(i,q_index) *
1.0 *
- fe_values.JxW(q_point));
+ fe_values.JxW(q_index));
}
// Finally, transfer the contributions from @p cell_matrix and