Re: [deal.II] InternalError in VectorTools::compute_nonzero_normal_flux_constraints

[Lars Corbijn van Willenswaard]

On Tue, Feb 21, 2017 at 11:28 PM, Wolfgang Bangerth 
wrote:

> On 02/21/2017 07:54 PM, Lars Corbijn van Willenswaard wrote:
>
>> I will try to do so somewhere this week, but it depends a bit on how much
>> free
>> evening time I have this week. Furthermore, I don't know what the
>> 'fall-out'
>> is of having to initialize the generalized support points, ass far as I
>> see
>> now I will have to implement and learn about the interpolation functions
>> of
>> FiniteElement.
>>
>
> I don't think anything particularly complicated is necessary -- you may
> just have to copy the support points from the base elements (or generalized
> points if the regular ones are not available) to the FESystem's ones. The
> interleaving will work just like for all of the other information being
> copied together.
>
> I have implemented building the generalized support points, which is
indeed not hard. I have not had enough time to look into the interpolation
functions. Those seem to require some more work to adjust the input values
before passing them to the base element, for which I have to study a bit on
how they work.


>
> For the change to VectorTools, I can do it. But I don't have an idea if the
>> code following it does assume in any way that it has normal support
>> points and
>> not generalized ones. If you (or somebody else with more experience) can
>> tell
>> me it doesn't I will also make a patch for it.
>>
>
> It's a good question whether that will work. I have to admit that I don't
> know -- it's quite possible that the function is only meant to be used for
> primitive elements. But the patch is definitely not going to make anything
> worse. Give it a try (in a separate pull request).
>

I als ochanged this, but it does not work. The problem is, that it trips
over 'face_system_to_component_index', which does not work for non
primitive elements (like the RT-elements I'm using). That function is used
to determine which sets of dofs form components at the support points.

Furthermore, even if this could be fixed for non primitive elements, I have
some doubts that the rest of the function will work with generalized face
support points and non primitive elements. Though I have to say that my
experience with either of them is practically zero. The algorithm that the
function follows is, if I read it correctly, to find at each support point
dim dofs, which form the dim components of the vector at that support
point. After handling the sharing of points/dofs between cells, it enforces
that the found dofs should match the normal constraint. Now I think that
for non primitive elements it will probably not find the dim independent
dofs which share the support point. While for generalized support points
the assumption that dofs can be coupled to a single support point and then
be constraint seems incorrect. Though again, this is speculation due to my
limited knowledge with both deal ii and non primitive elements.

Please let me know if I can do something with it, like putting a warning
above the relevant function. For my own project I think I will,
unfortunately, hate to use some alternative (hacky) ways to enforce the
constraint, due to time constraints.

Cheers,

Lars


>
> Cheers
>
>  W.
>
> --
> 
> Wolfgang Bangerth  email: bange...@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/fo
> rum/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/to
> pic/dealii/jFjGBuEhdXQ/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to output a single scalar in a parallel code

[Jean-Paul Pelteret]

Dear Pasha,

 I must correct my answer by replacing "residual vector" by "internal force 
> vector"


Ok, great. So this is a slightly different situation. I was going to say 
that I think that we're cross-talking here, and you saved me from sketching 
out a whole bunch of equations to clarify things :-) 

Yes, we can use internal forces for all situations. 


Hmm... I still don't think so. What about if you have body forces? Then the 
internal stresses balance both the externally applied loads plus the volume 
forces. Do you agree with this? However, one can of course compute the 
internal stress plus the body force loading... 

But then consider a cube with a non-uniform pressure on more than one side. 
How do I compute from this global vector the average or total load applied 
to any *individual* surface? You can't, because the values in the 
degrees-of-freedom shared by any two sides of the geometry contain load 
contributions from two faces. For reaction forces at constrained surfaces - 
no problem of course!

Also, I think that one should consider which vectors one typically 
constructs during assembly. The most efficient way to do assembly (but I 
recognise that this is not universally followed - not even by me) is to 
compute the constrained RHS vector and tangent matrix in one shot during 
assembly. That means that you never need to actually construct the 
individual internal, external and body force vectors in order to solve your 
linear system because, well, they get bundled into the RHS and then further 
manipulated when imposing constraints.

Dependent on what it is you want to achieve, I do concur that what you 
suggest may be perfectly plausible. But in general... I'm still 
unconvinced. To be introspective, what I suggest (and what Hamed has done 
here) might be overkill for plain elasticity alone since you can always 
recompute the traction vector itself as it is you that specified the 
loading at each time-/load-step. In coupled problems, however, external 
tractions are not always so easily defined and its actually sometimes 
easier to simply do this stress-type integration procedure.

However, can we calculate internal force by integrating over boundary 
> surfaces when material model is history dependent and has internal state 
> variable such as plasticity


Yes, you can - its just not very convenient. You could do this in one of 
two ways:

1. Extrapolate the internal variables to the faces and recompute the 
stresses with these values. Its clear that this might be problematic 
because one has to make an assumption about whether there is plastic yield 
at the extrapolated points (as can be attested to by the note in the 
ContinuousQuadratureDataTransfer 

 
class - we'd like to extend this at some point and we've already worked out 
a couple of options that may work for plasticity but simply haven't tried 
to implement them yet), so I would recommend option 2...
2. Simply create and track additional quadrature point data (with internal 
data) at the faces at which you'll be doing this post-processing. That way 
you can update the internal variables as necessary and they'll always be in 
sync with the solution field. You can then compute stresses etc. at these 
points whenever necessary. Its annoying if the local problem is expensive 
to compute, but the number of face QP data that you'd create is likely far 
less than those required for assembly etc.

Best,
Jean-Paul

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Access specific element within a distributed triangulation

['Seyed Ali Mohseni' via deal.II User Group]

Thank you everyone. I understand your suggestions. From what I learned now 
in deal.II is, that each cell is owned by a specific processor which means 
I cannot access the information within a locally owned cell. That's why I 
have to loop over all cells and work with distributed data without using 
locally owned entities. 

Maybe I describe my aim more properly. I like to find the maximum value 
within a MPI vector such as our locally_relevant_solution in step-40. 
Hence, I have to loop over all cells and all nodes to find where the 
geometrical position and value of the displacement lies. How is this 
possible within a parallely distributed environment?

My approach using local_dof_indices is somehow impossible. It only works 
for 1 core according to the below code:

std::vector displacement_norms, loaded_nodes, loaded_elements;
double max_displacement;
int max_loaded_node, max_loaded_element;

typename DoFHandler::active_cell_iterator cell = 
dof_handler.begin_active(), endc = dof_handler.end();
std::vector local_dof_indices(fe.dofs_per_cell);
unsigned int cell_number = 0;

for (; cell != endc; ++cell, ++cell_number)
{
 pcout << "CELL ID: " << cell_number << std::endl; // this works 
surprisingly.

 cell->get_dof_indices(local_dof_indices);// Output or accessing 
this won't work I assume.

 for (unsigned int i = 0, j = 1; i < fe.dofs_per_cell; i += 2, j += 2)
 {
  double displacement_x = locally_relevant_solution(local_dof_indices[i]);
  double displacement_y = locally_relevant_solution(local_dof_indices[j]);

  double displacement_norm = sqrt(displacement_x * displacement_x + 
displacement_y * displacement_y);

  loaded_nodes.push_back(cell->vertex_index(v));
  loaded_elements.push_back(cell_number);
  displacement_norms.push_back(displacement_norm);
 }
}

max_displacement = *max_element(displacement_norms.begin(), 
displacement_norms.end());

double max_index = distance(displacement_norms.begin(), max_element(
displacement_norms.begin(), displacement_norms.end()));

// Node and element number of maximal displacements
max_loaded_node = loaded_nodes[max_index];
max_loaded_element = loaded_elements[max_index];


Hopefully someone has an idea to help me out. Thank you :)


Kind regards,
S. A. Mohseni

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Access specific element within a distributed triangulation

['Seyed Ali Mohseni' via deal.II User Group]

Thank you everyone. I understand your suggestions. From what I learned now 
in deal.II is, that each cell is owned by a specific processor which means 
I cannot access the information within a locally owned cell. That's why I 
have to loop over all cells and work with distributed data without using 
locally owned entities. 

Maybe I describe my aim more properly. I like to find the maximum value 
within a MPI vector such as our locally_relevant_solution in step-40. 
Hence, I have to loop over all cells and all nodes to find where the 
geometrical position and value of the displacement lies. How is this 
possible within a parallely distributed environment?

My approach using local_dof_indices is somehow impossible:

for (; cell != endc; ++cell, ++cell_number)
{
 pcout << "CELL ID: " << cell_number << std::endl; // this works 
surprisingly.

 cell->get_dof_indices(local_dof_indices);// Output or accessing 
this won't work I assume.

 for (unsigned int i = 0, j = 1; i < fe.dofs_per_cell; i += 2, j += 2)
 {
  double displacement_x = locally_relevant_solution(local_dof_indices[i]);
  double displacement_y = locally_relevant_solution(local_dof_indices[j]);

double config_forces_norm = sqrt(config_forces_x * config_forces_x + 
config_forces_y * config_forces_y);
 }
}


-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to output a single scalar in a parallel code

[amir . pasha . beh]

Dear jean-Paul

Yes, we can use internal forces for all situations. I must correct my 
answer by replacing "residual vector" by "internal force vector". However, 
can we calculate internal force by integrating over boundary surfaces when 
material model is history dependent and has internal state variable such as 
plasticity? I think the answer is no, so I have suggested to use "internal 
force vector" in the first answer instead of integration over the boundary 
surfaces.


Thank you,
Pasha

On Sunday, February 26, 2017 at 5:06:18 PM UTC+3:30, Jean-Paul Pelteret 
wrote:
>
> Dear Amir,
>  
>
>> Why we must calculate reaction forces at a boundary which has external 
>> force? 
>>
>
> You certainly don't have to. In my example I used this as a demonstration 
> that the forces computed from the internal stresses matched the applied 
> (uniformly distributed) load. If you believe this then you could happily 
> apply the same procedure to any surface undergoing displacement control, 
> knowing that it will produce the desired result.
>  
>
>> But we usually calculate reaction forces at BCs with prescribed 
>> displacement which it is better to calculate them from residual vector 
>> (which equal to [-1.0 * internal forces] at these BCs). 
>>
> As the body is in equilibrium, they are balanced by the stresses generated 
>> by the material = internal forces
>>
>
> Ok, I see now where you're coming from. This might be true, but only for 
> constrained degrees-of-freedom for quasi-static problems. But in his 
> initial post Hamed wasn't specifically referring to reaction forces (unless 
> I misunderstood something), but rather the resultant force on an 
> unspecified boundary surface. Furthermore, for time-dependent problems you 
> may have additional contributions to these components of the residual due 
> to the imposition of constraints on the velocity and acceleration field (if 
> these constraints are non-homogeneous).
>  
>
>> Please, did you test the attached code for a plastic material model?
>>
>
> No I haven't. Do you anticipate the result will be any different when 
> using a plasticity model? If so, how?
>
> Regards,
> Jean-Paul
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to output a single scalar in a parallel code

[Jean-Paul Pelteret]

Dear Amir (and Hamed),
 

> Integration stress on external faces is correct when the material 
> constitutive model is linear (for example elastic material model in which 
> behavior is not history dependent). It is better to extract residual forces 
> at the boundary from the system residual vector. 
>

I disagree and consider both statements that you've made to be incorrect. 
For starters, in a nonlinear problem one uses a nonlinear solution scheme 
to drive the residual (a consequence of equilibrium imbalance) down to 
zero. So not only does the residual not represent any measure of external 
force, but also one expects this vector to be zero in for a balanced system!

In general (for a linear/nonlinear hyperelastic/dissipative material), to 
compute the total force acting on the surface of a domain one should 
integrate the tractions acting on that surface. As the body is in 
equilibrium, they are balanced by the stresses generated by the material. 
Using Cauchy's theorum 

 
one can re-express this traction in terms of the internal stresses and 
compute the forces developed on any arbitrary cut-plane in the material 
domain.

To demonstrate, I've attached an example for one-field elasticity that 
shows the computation of the reaction forces on either side of a cube under 
uniaxial stress conditions. All of the computations are done in the 
post-processing step (Solid::output_results). Since this cube has flat 
faces, computing the nodal reaction forces (which, by the way, are 
outputted so that you can visualise them as well) and then summing them is 
equivalent to directly integrating the contributions over the surface.

I hope that this helps provide you with some clarification on the point 
that you made.

Regards,
Jean-Paul

-- 
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.
For more options, visit https://groups.google.com/d/optout.
/* -
 * $Id: step-44.cc 30526 2013-08-29 20:06:27Z felix.gruber $
 *
 * Copyright (C) 2010 - 2013 by the deal.II authors and
 *  & Jean-Paul Pelteret and Andrew McBride
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -
 */

/*
 * Authors: Jean-Paul Pelteret, University of Cape Town,
 *  Andrew McBride, University of Erlangen-Nuremberg, 2010
 */


// We start by including all the necessary deal.II header files and some C++
// related ones. They have been discussed in detail in previous tutorial
// programs, so you need only refer to past tutorials for details.
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 


// We then stick everything that relates to this tutorial program into a
// namespace of its own, and import all the deal.II function and class names
// into it:
namespace Step44
{
  using namespace dealii;

// @sect3{Run-time parameters}
//
// There are several parameters that can be set in the code so we set up a
// ParameterHandler object to read in the choices at run-time.
  namespace Parameters
  {
// @sect4{Finite Element system}

// As mentioned in the introduction, a different order interpolation should be
// used for the displacement $\mathbf{u}$ than for the pressure
// $\widetilde{p}$ and the dilatation $\widetilde{J}$.  Choosing
// $\widetilde{p}$ and $\widetilde{J}$ as discontinuous (constant) functions
// at the element level leads to the mean-dilatation method. The discontinuous
// approximation allows $\widetilde{p}$ and $\widetilde{J}$ to be condensed
// out and a classical displacement based method is recovered.  Here we
// specify the polynomial order used to approximate the solution.  The
// quadrature order should be adjusted accordingly.
struct FESystem
{
  unsigned int poly_degree;
  unsigned int quad_order;

  static void