You are creating a new instance of FEValues at each quadrature point. This
is a very expensive operation, since there all the shape functions are
evaluated. Try to reuse that by passing it to the function as a parameter.
Hope that helps!
Peter
On Sunday, 8 January 2023 at 09:58:41 UTC+1 ce21...@smail.iitm.ac.in wrote:
> Following is my H_plus function:
>
> float PhaseField::H_plus(Vector<double> solution_elastic
> , const auto cell,const unsigned int q_point)
> {
> QGauss<2> quadrature_formula_damage(fe_damage.degree + 1);
>
> FEValues<2> fe_values_damage(fe_damage,
> quadrature_formula_damage,
> update_gradients |
> update_quadrature_points );
>
> fe_values_damage.reinit(cell);
>
> int node = 0;
>
> /* Initialising all strains as zero */
> float e_xx = 0.000;
> float e_yy = 0.000;
> float e_xy = 0.000;
>
> /*calculating strains*/
> for (const auto vertex : cell->vertex_indices())
> {
> int a = (cell->vertex_dof_index(vertex, 0));
> e_xx = e_xx +
> solution_elastic[a*2]*fe_values_damage.shape_grad(node, q_point)[0];
> e_yy = e_yy +
> solution_elastic[a*2+1]*fe_values_damage.shape_grad(node, q_point)[1];
> e_xy = e_xy
> +0.5*(solution_elastic[a*2]*fe_values_damage.shape_grad(node, q_point)[1]
>
> +solution_elastic[a*2+1]*fe_values_damage.shape_grad(node,
> q_point)[0]);
> node = node +1;
> }
>
> const auto &x_q = fe_values_damage.quadrature_point(q_point);
> float H_plus_val;
>
> /*This is the actual quantity I want to evaluate for each quadrature
> point*/
> H_plus_val = 0.5*lambda(x_q)*(pow((e_xx+e_yy),2))
> +
> mu(x_q)*((pow(e_xx,2))+2*(pow(e_xy,2)) + (pow(e_yy,2)));
>
> return H_plus_val;
> }
>
> H_plus function is called in assemble_damage for each quadrature point
>
> I am also attaching some lines of code from assemble_damage where H_plus
> is being called
>
> for (const auto &cell : dof_handler_damage.active_cell_iterators())
> {
>
> fe_values_damage.reinit(cell);
> cell_matrix_damage = 0;
> cell_rhs_damage = 0;
>
>
> float gc = 0.000027; //energy release rate
> float l = 0.015;
> float H;
>
> for (const unsigned int q_index :
> fe_values_damage.quadrature_point_indices())
> { float H_call = H_plus(solution_elastic,cell,q_index);
>
>
> if (H_call > H_vector[4*cell_number + q_index])
> {
> H = H_call;
> }
> else
> {
> H = H_vector[4*cell_number + q_index];
> }
> H_vector_new[4*cell_number + q_index] = H;
> for (const unsigned int i : fe_values_damage.dof_indices())
> {
>
>
> for (const unsigned int j : fe_values_damage.dof_indices())
> {
>
> const auto &x_q = fe_values_damage.quadrature_point(q_index);
>
>
> cell_matrix_damage(i, j) +=
> // contribution to stiffness from -laplace u term
>
> Conductivity_damage(x_q)*fe_values_damage.shape_grad(i, q_index) * //
> kappa*grad phi_i(x_q)
> fe_values_damage.shape_grad(j, q_index) * // grad
> phi_j(x_q)
> fe_values_damage.JxW(q_index) // dx
> +
> // Contribution to stiffness from u term
>
>
> ((1+(2*l*H)/gc)*(1/pow(l,2))*fe_values_damage.shape_value(i, q_index) *
> // phi_i(x_q)
> fe_values_damage.shape_value(j, q_index) * //
> phi_j(x_q)
> fe_values_damage.JxW(q_index)); //
> dx
> }
>
> cell_rhs_damage(i) += (fe_values_damage.shape_value(i, q_index) *
> // phi_i(x_q)
> (2/(l*gc))* H*
> fe_values_damage.JxW(q_index)); // dx
> }
> }
>
> /*Adding the local k and local f to global k and global f*/
> cell->get_dof_indices(local_dof_indices);
> for (const unsigned int i : fe_values_damage.dof_indices())
> {
> for (const unsigned int j : fe_values_damage.dof_indices())
> system_matrix_damage.add(local_dof_indices[i],
> local_dof_indices[j],
> cell_matrix_damage(i, j));
>
> system_rhs_damage(local_dof_indices[i]) += cell_rhs_damage(i);
> }
> cell_number = cell_number + 1;
>
> }
>
> Thanks and regards
> Wasim
>
> On Saturday, January 7, 2023 at 11:59:49 PM UTC+5:30 blais...@gmail.com
> wrote:
>
>> There might be many things that can be done to improve the speed of this
>> function. You can ask yourselve the following question as guidance:
>> - Does the function allocate memory?
>> - Could it be inlined?
>> - Are you calling the function inside the DOF loops or inside the
>> quadrature loop?
>>
>> Then I would time the function to measure if this is actually the real
>> culprit or if it could be something else.
>> If you copy/paste the content of your assembly code and the function, I
>> would be glad to give it a look (and I am sure others here will help you
>> too).
>>
>>
>> On Saturday, January 7, 2023 at 12:02:13 a.m. UTC-5
>> ce21...@smail.iitm.ac.in wrote:
>>
>>> Sorry for the confusion. I think I made a mistake while writing the
>>> first email.
>>>
>>> H_plus is being called in Assemble_damage and not assemble_elastic. It
>>> uses elastic solution, cell and gauss point to evaluate strain at a gauss
>>> point. Then some quantity is evaluated based on the strain.
>>>
>>> Similarly I have another function damage_gauss which is being called in
>>> assemble_elastic that evaluates damage at a gauss point using the damage
>>> solution, cell and gauss point.
>>>
>>> Wasim Niyaz
>>> Research scholar
>>> CE Dept.
>>> IITM
>>>
>>> On Sat, 7 Jan, 2023, 10:15 am Wasim Niyaz Munshi ce21d400, <
>>> ce21...@smail.iitm.ac.in> wrote:
>>>
>>>> I use it to evaluate strain at Gauss points. Then, i evaluate some
>>>> quantity which is a function of this strain.
>>>>
>>>> Wasim Niyaz
>>>> Research scholar
>>>> CE Dept.
>>>> IITM
>>>>
>>>> On Sat, 7 Jan, 2023, 3:09 am Wolfgang Bangerth, <bang...@colostate.edu>
>>>> wrote:
>>>>
>>>>> On 1/6/23 13:53, Wasim Niyaz Munshi ce21d400 wrote:
>>>>> > I am using 65536 elements. For step-8 the assembly takes very less
>>>>> time
>>>>> > (around 0.15second) while for my assemble_elastic, it takes around 5
>>>>> seconds.
>>>>> > The only difference between my assemble_elastic function and the
>>>>> assemble
>>>>> > function of step-8 is that for each Gauss point, I additionally
>>>>> call a
>>>>> > function(H_plus) that takes the laplace solution, the current cell
>>>>> and Gauss
>>>>> > point as input and evaluates some quantity using this information.
>>>>> > The H_plus function is called 4*65536 times but the function is very
>>>>> simple.
>>>>> > My doubt is whether such a huge increase in cost (from 0.15 sec to 5
>>>>> sec) is
>>>>> > expected for this problem or is there something that I am doing that
>>>>> is
>>>>> > increasing the cost so much.
>>>>>
>>>>> Wasim, the question is what you do with "the laplace solution, the
>>>>> current
>>>>> cell and Gauss point". If you show us what H_plus does, we may be able
>>>>> to advise.
>>>>>
>>>>> Best
>>>>> W.
>>>>>
>>>>> --
>>>>>
>>>>> ------------------------------------------------------------------------
>>>>> Wolfgang Bangerth email:
>>>>> bang...@colostate.edu
>>>>> www:
>>>>> http://www.math.colostate.edu/~bangerth/
>>>>>
>>>>>
>>>>> --
>>>>> The deal.II project is located at http://www.dealii.org/
>>>>> For mailing list/forum options, see
>>>>> https://groups.google.com/d/forum/dealii?hl=en
>>>>> ---
>>>>> You received this message because you are subscribed to a topic in the
>>>>> Google Groups "deal.II User Group" group.
>>>>> To unsubscribe from this topic, visit
>>>>> https://groups.google.com/d/topic/dealii/-6ndTW_k5fQ/unsubscribe.
>>>>> To unsubscribe from this group and all its topics, send an email to
>>>>> dealii+un...@googlegroups.com.
>>>>> To view this discussion on the web visit
>>>>> https://groups.google.com/d/msgid/dealii/895a079c-2b85-b14f-94ee-b4b78336884d%40colostate.edu
>>>>> .
>>>>>
>>>>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/bc38adc5-563b-476a-82c4-255246af25a5n%40googlegroups.com.
