________________________________ From: Andrew McBride <[email protected]> To: Pietro Maximoff <[email protected]> Sent: Tuesday, September 29, 2009 9:56:42 AM Subject: Fwd: [deal.II] 1D Boundary-Value Problem see Mike's comments below. Begin forwarded message: From: Michael Rapson <[email protected]> >Date: 29 September 2009 10:25:16 AM >To: Andrew McBride <[email protected]> >Subject: Re: [deal.II] 1D Boundary-Value Problem > > >Hi there Pietro and Andrew, > >Two points about the code given below: the loop over vertices seems to >be unnecessary, and due to c style numbering fe.dofs_per_cell will >give you a value too large by one for cell_rhs. My two cents: > > for (; cell!=endc; ++cell) > { > fe_values.reinit (cell); > cell_matrix = 0; > cell_rhs = 0; > > for (unsigned int q_point=0; q_point<n_q_points; ++q_point) > > { > 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)) * fe_values.quadrature_point(q_point)[0]; > > 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)); > } > } > > > if (cell->at_boundary(1) == true) > { > const double neumann_value = -0.5; > cell_rhs(1) += (neumann_value); // see comment below on >why 1 is used > } > } // end of loop on cells > > // Copy local matrix into global matrix > > >In the code above I use the original test for a cell being at a >boundary, in the loop over cells which will iterate over all cells. If >we are at the end cell the if statement executes. There are further >complexities in choosing the correct dof to assign the neumann_value >to because the ordering of dofs in the FE is defined by the library >internals only and should not be relied on. For the FE_Q class the >documentation specifies that LHS dof gets 0 and RHS dof gets 1, the >remaining dofs are assigned to any internal points. Hence 1 seems to >be appropriate for FE_Q in this case but writing code that depends on >this ordering is not advised. I cannot think of the appropriate method >to use to avoid it in this case however. > >Cheers, >Michael > >Hi Michael and Andrew Fantastic. Works like a charm. Thanks . Pietro
_______________________________________________ dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
