Hi, The integral is an 'average' weighted by the size of the cell. (Think about the integral of a constant, say t, from 0 to 1 vs from 0 to 0.5. The second 'cell' is half the size of the first and the integral is half the value.) This is desirable behavior in some cases, for example if you were summing each integral over cells over the entire domain, but not in other cases, e.g. if you want to compare the integral on one cell to another on a cell of a different size.
If you want the size of the cell to play a role you can leave it as is, dividing by the cell volume will give you the second type of behavior. Regards, Michael On Mon, May 17, 2010 at 11:56 AM, Isabel Gil <[email protected]> wrote: > Hi! > > Thanks for your answer. > > Andrew McBride escribió: > > Hi Isabel > > I'm confused by a couple of things. > > You want to determine the trace of the Laplacian of Phi at the quadrature > points and then extrapolate the calculated values back to the nodes? Why are > you multiplying the value at the quadrature point by fe_values.JxW(q_point)? > This implies that you are doing integration based on nodal values as opposed > to extrapolation. What you want to do is just extrapolate the quadrature > point data back to the nodes without integrating. The function you need to > do this FETools::compute_projection_from_quadrature_points_matrix > The projection is on an element so the result will be globally > discontinuous. > > > What I want to calculate is the Laplacian of Phi over the nodes of my mesh. > > When I use the loop "for (unsigned int > q_point=0;q_point<n_q_points;q_point++)", I always multiply by > "fe_values.JxW(q_point)" to get the transformation from reference to real > cell. > > In calculating the integral why are you using the > fe_face_values.get_function_values (rho,values)? I would extrapolate the > values at the nodes to the quadrature points and then integrate them using a > sum over the quadrature points and multiplying each quadrature point > contribution by fe_values.JxW(q_point). There may even be a deal function > that does this? > > > Sorry, I made a mistake :-[ . What I wanted to use (and indeed I use) was > "fe_values.get_function_values (rho,values);" > But, once I have the "integral" value, should I divide it by the volume to > get the average or does the "integral" value itself contain the avarage? > > Thanks again. > Best. > Isa > > Cheers > Andrew > > > > Am 17 May 2010 um 10:58 AM schrieb Isabel Gil: > > > > Hi! > > I have the value of the scalar phi(x,y,z) and I would like to calculate > rho(x,y,z)=- const*div(grad(phi)), so I do: > . > . > . > for (; cell_p!=endc; ++cell_p) > { > cell_rhs=0.0; > fe_values.reinit (cell_p); > fe_values.get_function_2nd_ > derivatives (phi, values); > for (unsigned int i=0; i<dofs_per_cell; ++i) > > cell_rhs(i)-=const*fe_values.shape_value(i,q_point)* > > (values[q_point][0][0]+values[q_point][1][1]+values[q_point][2][2])*fe_values.JxW(q_point); > > cell_p-> get_dof_indices(local_dof_indices); > for (unsigned int i=0;i<dofs_per_cell;++i) > rho(local_dof_indices[i]) += cell_rhs(i); } > > Secondly, I would like to calculate the average rho, so I do: > > double integral=0.0; > > for (; cell_p!=endc; ++cell_p) > { > fe_values.reinit (cell_p); > fe_face_values.get_function_values (rho,values); > > for (unsigned int q_point=0;q_point<n_q_points;q_point++) > integral += values[q] * fe_face_values.JxW(q); > } > > Are both pieces of code all right? > > Does the value of integral correspond to {(1/Volume)*integral (rho)} or only > to {integral(rho)}? > > Thanks in advance. > Best > Isa > > _______________________________________________ > dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii > > > > > > _______________________________________________ > dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii > > _______________________________________________ dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
