Re: [deal.II] Mesh-induced elastic anisotropy and distorting the quad. points as a way to palliate it

[David F]

Hi Sebastian,

first of all, sorry for my late reply. Thank you very much for comment, it 
certainly raises some very interesting points. I think the only thing which 
is still left to be explained is the fact that a triangular grid, which 
yields the right isotropic result, becomes anisotropic just by rearranging 
the elastic properties to form squares-like clusters. I think this case 
does not correspond to any of your 4 points since I do this without h->inf. 
In fact, each of those clusters are just 2 FE.

Best,
David.


On Saturday, 11 July 2020 06:19:50 UTC+2, Sebastian Stark wrote:
>
> Hi David,
>
> I'm neither an expert, nor do I know the literature well, but looking on 
> your pictures, I think, the situations you are studying are geometrically 
> anisotropic. Just plot the distribution of angles the faces of your 
> inhomogenities make with the x-axis. For the quad-case, you'll get two 
> discrete peaks at 0 and 90 degree. for the triangular case, you get 0, 45 
> and 90 degree. So, from this, the results do not seem surprising to me 
> (just consider the extreme case of cracks - if you have them only at 0 and 
> 90 degree oriented, this is unlikely to be isotropic). The fact that you 
> have randomly assigned elastic properties won't help to fix that.
>
> A few examples (in 2d):
>
> (1) equidistant circular inclusions in a matrix (matrix and inclusions two 
> different isotropic linearly elastic materials) -> this should be isotropic 
> if mesh size h->0
>
> (2) equidistant square inclusions in a matrix, all aligned with x-axis 
> (matrix and inclusions two different isotropic linearly elastic materials) 
> -> probably anisotropic if mesh size h->0
>
> (3) equidistant square inclusions in a matrix, random orientation of 
> inclusions (matrix and inclusions two different isotropic linearly elastic 
> materials) -> should be isotropic if mesh size h->0 and number of 
> orientations->infinity
>
> (4) equidistant square inclusions in a matrix, random orientation of 
> inclusions (random isotropic linearly elastic materials) -> should be 
> isotropic if mesh size h->0 and number of orientations->infinity
>
> Also to consider: In 3d, there are 21 elastic constants for a linearly 
> elastic material. In a mathematical 2d scenario, it should be 6. This 
> suggests that, in example (3), it is not strictly necessary to have random 
> orientation. Rather, a few (equally spaced) discrete orientations might be 
> good enough. If that's the case, how many does one need? I'm betting on 6, 
> not sure though. Alternatively, one could replace the square inclusions in 
> (2) by regular polygons and ask how many vertices the polygon needs for 
> isotropy. Again, I'm betting on 6.
>
> Related: Are there crystal structures with such a high degree of symmetry, 
> that they are elastically isotropic? For dielectric properties, a cubic 
> crystal is good enough already. But the dielectric tensor is rank 2 and the 
> elastic one rank 4. So you'll need much more symmetry in the crystal; and 
> considering that a crystal can have 6-fold rotational symmetry at most that 
> might be impossible.
>
> What I'm just noticing: Hexagonal crystals are elastically isotropic 
> perpendicular to the hexagonal axis. So, my bet on 6 might be good. And it 
> might explain your observation that the triangular elements are relatively 
> isotropic (though maybe not perfectly).
>
> I hope that gives you some input. If you have definitive answers to any of 
> the questions, I'm curious.
>
> Regards,
> Sebastian
>
>
> Am 11.07.20 um 00:12 schrieb David F:
>
> I have made a somewhat extensive study on his issue and prepared some 
> plots that will hopefully answer your questions, and also includes Bruno's 
> suggestion about distorting the mesh. The basic setup is: I sheared the 
> mesh along different orientations (see x-axis on the plots) and measured 
> the shear modulus (y-axis). I have repeated the random process of setting 
> the elastic properties many times to have good statistics (see errorbars on 
> the plots). Each element has an isotropic stiffness tensor with a Poisson 
> ratio of 1/3 and a shear modulus which is exponentially distributed with an 
> average of 10. I use linear shape functions unless otherwise stated. If the 
> picture are not big enough, you can find them in the links beneath them.
>
>
> *1) I change the resolution. *By this I don't mean just a mesh with a 
> bigger number of elements, but importantly each inhomogeneities is 
> represented by a bigger number of elements. Therefore, we solve problems 
> with exactly the same physical domain but with different resolution. In the 
> legend, n me

Re: [deal.II] Mesh-induced elastic anisotropy and distorting the quad. points as a way to palliate it

[David F]

I have a 2D system for which I create the stiffness tensor of an isotropic 
material, but for each finite element I create it with a different shear 
modulus. The shear modulus is random for each element (I use an exponential 
distribution, but any distribution leads to the same behavior as long as 
the std is high), with no structure such as layers or anything else. In 
this case, the system should clearly be macroscopically isotropic (up to 
statistical fluctuations due to the random properties) for symmetry reasons.

On Thursday, 9 July 2020 05:03:45 UTC+2, Wolfgang Bangerth wrote:
>
> On 7/2/20 10:06 PM, David F wrote: 
> > 
> > *_Q2_:* why the system behaves as anisotropic if its local inhomogeneous 
> > elastic properties are isotropic? If you have any comment or suggestion 
> about 
> > the problem of mesh-induced elastic anistropy in FEM, I would like to 
> know it. 
>
> I don't know how exactly you choose your coefficient, but if you alternate 
> layers of isotropic materials, then you get an anisotropic material. Think 
> about layering styrofoam plates with steel plates -- the resulting stack 
> of 
> layers is essentially incompressible under loads from the side (because 
> the 
> steel plates provide the stiffness), but is quite compressible if you load 
> it 
> from top and bottom (because the styrofoam layers will simply collapse). 
>
> 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 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/ef939354-6b36-4b8d-aaa8-f15e5d6aca3bo%40googlegroups.com.

[deal.II] Mesh-induced elastic anisotropy and distorting the quad. points as a way to palliate it

[David F]

Hello everyone,

I'm trying to solve a 2D solid mechanics homogenization problem, in which I 
have element-wise constant elastic properties, which are inhomogeneous and 
isotropic from element to element (i.e., I am assembling the system using 
the same 4-rank stiffness tensor for all the quadrature points of a certain 
element, but that tensor is different for each element). For this system, I 
would like to compute its effective elastic properties, which I do by 
loading the system under several different loading conditions as is done in 
standard homogenization approaches. The system should behave as an 
isotropic solid. However, I observe significant anisotropy (and clearly not 
due to random fluctuations that might arise because the element-to-element 
inhomogeneous properties are randomly distributed). I attribute this to a 
mesh dependency of the solution, since I have solved the same problem with 
a unstructured triangular mesh with another FEM package and I don't observe 
this issue. I believe the structured quadrilateral mesh induces some 
artificial elastic anisotropy, which is not there in the case of the 
unstructured triangular mesh due to its topological disorder.

I've thought of a way that might palliate this issue, which is to set 
different elastic properties at the quadrature points themselves (i.e., the 
properties are no longer element-wise constant). This seems to work to some 
extent since the system becomes less anisotropic, however it is not good 
enough.


*Q1:* is there a preferred way in dealII in which I could randomly distort 
a bit the location of the quadrature points? I think this extra distortion 
might help get rid of the mesh artifacts. Is is possible to do it with the 
in-built Lagrange linear FE or another type of FE is more suitable within 
dealII for this task? Basically I have no idea where to start from to do 
something like this, so any suggestion is welcome.

*Q2:* why the system behaves as anisotropic if its local inhomogeneous 
elastic properties are isotropic? If you have any comment or suggestion 
about the problem of mesh-induced elastic anistropy in FEM, I would like to 
know it.


Thanks in advance,
David.

-- 
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/a6e817f2-7244-4717-b1e2-f112e7ce6888o%40googlegroups.com.

Re: [deal.II] Using different random seed when distorting a triangulation

[David F]

Yes, it seems indeed easy. I'll try to learn how to send it through the 
official channel in github.

On Tuesday, 26 May 2020 05:40:57 UTC+2, Wolfgang Bangerth wrote:
>
> On 5/25/20 9:20 PM, David F wrote: 
> > 
> > I would like to apply many random distortions to a triangulation. 
> However, the 
> > function GridTools::distort_random seems to produce always the same 
> > distortion. Is there a way to use different seeds in order to generate 
> always 
> > different distortions? It seems that it's not possible through the 
> interface 
> > of GridTools::distort_random but I wonder if there is another way to do 
> it, 
> > since this seems to me a very basic feature for the users of this 
> function. 
>
> Interesting question! 
>
> The function doesn't currently allow that, it uses the same sequence of 
> random 
> numbers every time -- by design: This way, if you run your program twice 
> in a 
> row, it'll produce the same result. 
>
> But it wouldn't be very difficult to change that. The function is here: 
>
> https://github.com/dealii/dealii/blob/master/source/grid/grid_tools.cc#L1013-L1244
>  
> The random number generator is initialized in line 1099. It wouldn't be 
> very 
> difficult to add an argument to the function that acts as a random seed 
> and is 
> used in that line to seed the boost::random::mt19937 generator. Want to 
> give 
> this a try and implement such a change? We'd of course be very happy to 
> take a 
> patch! 
>
> Best 
>   Wolfgang 
>
>
> -- 
>  
> 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 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/1d6b5b3f-03b4-4de0-b186-f7ca9612f2a2%40googlegroups.com.

[deal.II] Using different random seed when distorting a triangulation

[David F]

Hello everyone.

I would like to apply many random distortions to a triangulation. However, 
the function GridTools::distort_random seems to produce always the same 
distortion. Is there a way to use different seeds in order to generate 
always different distortions? It seems that it's not possible through the 
interface of GridTools::distort_random but I wonder if there is another way 
to do it, since this seems to me a very basic feature for the users of this 
function.


Best,
David.

-- 
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/a08df296-b37f-40b4-8b4d-05f7126c22d8%40googlegroups.com.

[deal.II] Re: Periodicity with simultaneous displacement controlled BCs

[David F]

Thank you all, I found the way to do it. For future users trying to solve 
the same problem, I also found extremely helpful this paper:

*Effective properties of composite materials with periodic microstructure: 
a computational approach, J.C. Michel, H. Moulinec, P. Suquet, Comput. 
Methods Appl. Mech. Engrg. 172 (1999) 109-143*


David.

-- 
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/804705b2-e31e-4dec-aea8-f700e17d3140%40googlegroups.com.

[deal.II] Periodicity with simultaneous displacement controlled BCs

[David F]

Hello everyone, 

I'm trying to compute the effective elastic properties of a heterogeneous, 
linear and a bi-periodic system (i.e., left-right and top-bottom periodic 
displacement fields). To this system, I would like to apply a global 
shearing by prescribing the displacement field of the surfaces in the form 
of Dirichlet BCs. This seems slightly contradictory since a bi-periodic 
system doesn't have surfaces.

However, we can still think of the global shearing as an average surface 
displacement around which periodic fluctuations occur (the origin of such 
fluctuations is due to heterogeneous elastic properties). This is 
illustrated in the picture below.


[image: bitmap.png]

I would like to know which is the best way to do this in deal.II (I have 
tried with make_periodicity_constraints and interpolate_boundary_values, 
but the problem is that, as I explain before, we set apparently 
contradictory constraints).

Thank you in advance.


P.S.: suggestions of alternative procedures to achieve the same goal are 
welcome.

-- 
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/0c2893c7-e627-42ed-873f-38415a4c78ab%40googlegroups.com.

Re: [deal.II] Is parallel direct solver extremely slow?

[David F]

Dear Pai, 

I'm very interested in solving a problem with characteristics very similar 
to yours. Consequently, I run your modified code of step-17.cc for 30*30*30 
cells and for me it takes 6.43s with cg with -np 2 instead of 0.39s. Do you 
have any idea where this huge speed up migth come from? Is it maybe due 
some optimezed libraries that you are using? Right now I'm using the system 
default's one that I can find in the OS repositories. I would really 
appreciated if you could give me some hint about this, or which strategy 
you found more effective for solving many times the same elastic system 
with different rhs.

Best regardsm
David.

On Saturday, 1 September 2018 04:56:32 UTC+2, Pai Liu wrote:
>
> Hi Wolfgang,
>
> Thank you so much for all your detailed explanation. Now I have a general 
> idea of what all these things are and what I should do for my problem (a 
> multiple load case problem). I really appreciate your kind help.
>
> BlockJacobi builds an LU decomposition of that part of the matrix that is 
>> stored locally. So it's a really expensive preconditioner to build (which 
>> I 
>> gather you don't include in the time?) but then solves the problem in 
>> only a 
>> few iterations. If you want a fair comparison, you need to include the 
>> time to 
>> build the preconditioner. 
>
> However, I really like to figure out the problem of timing you mentioned. 
> Here I attach the step-17.cc I modified. 
> *I just added timing and modify the meshing code (to generate 30*30*30 
> cells), and nothing else.*
> I added timing to the member function solve() like the following and 
> change nothing else.
>
>   template 
>   unsigned int ElasticProblem::solve ()
>   {
> *TimerOutput::Scope t(computing_timer, "solve");*
>
> SolverControl   solver_control (solution.size(),
> 1e-8*system_rhs.l2_norm());
> PETScWrappers::SolverCG cg (solver_control,
> mpi_communicator);
>
> PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
>
> cg.solve (system_matrix, solution, system_rhs,
>   preconditioner);
>
> Vector localized_solution (solution);
>
> hanging_node_constraints.distribute (localized_solution);
>
> solution = localized_solution;
>
> return solver_control.last_step();
>   }
>
> *Thus I think the timing includes both the time to build the 
> preconditoner, as well as the time to solve Ax=b;*
>
> *And when I run this file with mpirun -np 2 ./step-17, it really just 
> takes about 0.4s to solve a 30*30*30 cells problem (with all the boundary 
> conditions unchanged from the original step-17 example).*
>
> Best,
> Pai
>

-- 
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.

Re: [deal.II] Obtaining the final position of the vertices

[David F]

Thank you very much. With those functions in DofTools I should be able to 
do what I need.

Best,
David.

On Tuesday, 19 March 2019 17:41:32 UTC+1, Wolfgang Bangerth wrote:
>
> On 3/19/19 9:57 AM, David F wrote: 
> > 
> > I am not sure how to answer your question. I'm using a very basic a 
> setup 
> > equivalent to step-8. Therefore, I have a solution vector with final 
> > displacements where each entry corresponds to the displacement of a dof. 
> My 
> > aim is to find the initial position of the vertices in the form of 
> > std::vector> and the final position (i.e., the initial 
> position + 
> > its displacement) in the same form. 
>
> I see -- you are using the pre-existing mesh and want to interpret the 
> solution as the *displacement* to the mesh. 
>
> If all you are interested in is to output the solution on the deformed 
> mesh, 
> you can call Data::build_patches() with the mapping argument, where you 
> use a 
> MappingQEulerian object initialized with your displacement vector. 
>
> If you want to have a representation in the program, you can use the 
> function 
> in namespace DoFTools that returns a vector of support point locations for 
> each DoF, and then query your solution vector for the corresponding 
> displacements. You might have to juggle with the x, y, and z displacements 
> a 
> bit to get this right, but the infrastructure is all there. 
>
> 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 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.

Re: [deal.II] Obtaining the final position of the vertices

[David F]

Dear prof. Bangerth, 

I am not sure how to answer your question. I'm using a very basic a setup 
equivalent to step-8. Therefore, I have a solution vector with final 
displacements where each entry corresponds to the displacement of a dof. My 
aim is to find the initial position of the vertices in the form of 
std::vector> and the final position (i.e., the initial position 
+ its displacement) in the same form.


Best,
David.

On Tuesday, 19 March 2019 16:17:58 UTC+1, Wolfgang Bangerth wrote:
>
> On 3/19/19 9:12 AM, David F wrote: 
> > 
> > I want to obtain the final position of the vertices (specifically, the 
> > vertices at the faces), i.e., the deformed configuration. I think that a 
> way 
> > of doing this is by creating a set of points and using fe_values to 
> > extrapolate the solution to those points. However, I don't like the idea 
> of 
> > extrapolating since the displacements are actually known in the vertices 
> > already. Could somebody give me some hint about what's the best way to 
> do this? 
>
> David -- how do you actually deform the mesh? Are you moving around 
> vertices? 
> In that case, you know where you move them to. Or are you using a 
> MappingEulerian? 
>
> 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 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] Obtaining the final position of the vertices

[David F]

Hi all,

I want to obtain the final position of the vertices (specifically, the 
vertices at the faces), i.e., the deformed configuration. I think that a 
way of doing this is by creating a set of points and using fe_values to 
extrapolate the solution to those points. However, I don't like the idea of 
extrapolating since the displacements are actually known in the vertices 
already. Could somebody give me some hint about what's the best way to do 
this?

Thanks in advance.

-- 
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] Calculating gradients at arbitrary points

[David F]

Hello, I'm trying to calculate the symmetric gradients of a solution vector 
at an arbitrary set of points. However, so far I have only found the way to 
do it at the quadrature points by means of 
*fe_values[fe_extractor].get_function_symmetric_gradient*. 

I would like to know what is the proper way to do it (calculate the 
symmetric gradients at an arbitrary set of points) in dealII at any other 
point (I mean using as much as possible dealII's methods), and if so, is 
there's a performance difference between calculating it at quadrature 
points and at another set of points (obviously, comparing set of same size)?

Thank you very much for your help, David.

-- 
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.

Re: [deal.II] Parasitic stress in corners with periodic boundary conditions

[David F]

Hi, sorry for this late reply. I missed it and realized just now. I found 
the problem, which is simply that if the whole system is fully periodic ( I 
mean, for example from plane to torus topology) then the 4 corners become 
effectively fixed (the other possibility would be rigid body motion). The 
fact that this 4 nodes become fixed, makes that when some eigenstrain is 
prescribed somewhere, unexpected stress appears close to the corners as the 
system can't fully relax there. I think this is something one was to live 
with. I didn't see it with other FEM packages because the solution suffered 
some post processing before I got it, which smoothed the result and somehow 
hide this.


On Monday, 2 June 2014 00:59:11 UTC+2, Wolfgang Bangerth wrote:
>
> On 05/28/2014 10:29 AM, David F wrote: 
> > Hello, my problem is the following: 
> > 
> > I prescribe an eigenstrain value in one element of the grid (i.e., an 
> inner 
> > element undergoes a certain transformation and this introduces strain 
> and 
> > stress in the rest of the grid), and everything works fine for normal 
> > boundaries and periodic boundaries. However, if one pays close 
> attention, for 
> > the periodic boundary solution in the corners of the grid one can see 
> some 
> > parasitic stress that shouldn't be there. 
> > 
> > I think the reason is the following: prescribing, let's say, a shear 
> > eigenstrain in one element, the four corners (in 2D) of the grid deform 
> in 
> > opositte directions if the boundary is not periodic, therefore when it 
> is 
> > periodic the four corners are in fact "the same", and these 4 opposite 
> > displacements would add up to 0, effectively fixing the "corner node". 
> This 
> > seems to me a natural problem of FEM and periodicity, so I have no clue 
> how to 
> > correct this with code. 
> > 
> > I have used other FEM packages for solving exactly the same problem I 
> > described, and somehow this is not happening, so there must be a way to 
> avoid 
> > this. I am not sure if it has to do with my implementation of the 
> periodic 
> > boundaries, or maybe the way in which deal.ii deals with periodicity is 
> the 
> > reason, or maybe I have to apply some kind of postprocessing correction 
> of 
> > which I am not aware yet. 
> > 
> > Does anyone know why could this be? 
>
> I don't know myself. Maybe someone else does. But I'd like to point out 
> that 
> using periodic boundary conditions is equivalent to looking at a problem 
> where 
> you have periodic array of sources. Is the element where you prescribe the 
> eigenstrain by chance close to a corner of your domain? If so, I would 
> expect 
> there to be some strain close to the opposite corner as well, simply 
> because 
> there is a source (outside your domain) close to that opposite corner. 
>
> It may be worthwhile sending a picture that shows what you are observing. 
>
>
> > P.S.: my periodic boundary implementation is a direct extension of 
> step-45 and 
> > everything but these small parasitic stress in the corners is pretty 
> accurate, 
> > so probably a problem in the code is not reason. Should I try with 
> > DoFTools::make_periodicity_constraints or do you the problem will be 
> exactly 
> > the same? 
>
> Hard to tell. It's probably worth a try. 
>
> Best 
>   W. 
> -- 
>  
> Wolfgang Bangerth   email:bang...@math.tamu.edu 
>  
>  www: http://www.math.tamu.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 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.

Re: [deal.II] matrix factorization reuse?

[David F]

Hi, I'm interested in this as well. Could you post or link to your 
solution? Thanks.


On Friday, 19 February 2016 17:21:22 UTC+1, Michael Harmon wrote:
>
> If you want to use Trilinos, I refactored their direct solver so that it 
> can be used the way you want, I will submit the patch, but if you want can 
> send you the changes as well.
>
>
> On Tuesday, December 1, 2015 at 4:05:29 AM UTC-6, Pascal Kraft wrote:
>>
>> Hello Wolfgang,
>>
>> does the same hold true for other direct solvers? I would like to use 
>> Mumps (PETSc) to get a LU-factorization and reuse that a couple of times. 
>> The PETSc wrappers only provide a constructor and solve. If I reuse solve 
>> withe the same matrix argument, will the factorization be done again or 
>> does it detect that the inverse is already present? The code for both the 
>> PETSc wrappers as well as the Trilinos Wrappers does not look as if that 
>> was supported.
>>
>> Greetings,
>> Pascal Kraft
>>
>> Am Sonntag, 4. Oktober 2015 20:14:30 UTC+2 schrieb bangerth:
>>>
>>>
>>> Zhen, 
>>>
>>> Jean-Paul already gave the correct answer, namely use a sparse direct 
>>> solver 
>>> and factorize only once. 
>>>
>>> > Therefore my following questions are: 
>>> > 
>>> > What is the limit of the size of matrix system can I do a lapack LU 
>>> > factorization on a single machine? 
>>>
>>> For a sparse direct solver, I think you can probably go to around 1M 
>>> unknowns. 
>>> If you use a LAPACK matrix, you will quickly be limited by memory. For 
>>> example, a 10,000 x 10,000 matrix already requires almost a GB of 
>>> memory. 
>>>
>>>
>>> > Are there any better way to get LU factorization in deal.ii? 
>>> > Can we incorporate "history experience" to achieve more efficient 
>>> solvers? 
>>>
>>> People have tried that and there is a substantial amount of literature 
>>> on the 
>>> topic of Krylov subspace recycling. I always think that it is a great 
>>> idea and 
>>> want to investigate using it, but it has not found its way into 
>>> mainstream 
>>> applications to the best of my knowledge. 
>>>
>>> Best 
>>>   W. 
>>>
>>> -- 
>>>  
>>> Wolfgang Bangerth   email:bang...@math.tamu.edu 
>>>  www: 
>>> http://www.math.tamu.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 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: PETScWrappers::MPI::Vector Vs parallel::distributed::Vector

[David F]

Hi, thanks for your answer. I have measure the time it takes for 
PETScWrappers::MPI::Vector, parallel::distributed::Vector< Number > and 
Vector< Number > to complete a very simple task consisting of accessing the 
elements of these vectors. Something like this (repeating this whole 
process 15 times for averaging results, and using very big vector sizes):

double a;
for(unsigned int i=0; i
> Hi,
>
> I guess it's more a question of preference. What I do is using the same 
> vector type as the matrix type: PETSc matrix -> PETSc vector, Trilinos 
> matrix -> Trilinos vector, matrix-free -> deal.II vector. deal.II vector 
> can use multithreading unlike PETSc vector but if you are using MPI, I 
> don't think that you will see a big difference.
>
> Best,
>
> Bruno
>
> On Thursday, August 25, 2016 at 5:30:31 PM UTC-4, David F wrote:
>>
>> Hello, I would like to know if among PETScWrappers::MPI::Vector and 
>> parallel::distributed::Vector< Number >, one of them is preferred over the 
>> other. They both seem to have a similar functionality and a similar 
>> interface. Although parallel::distributed::Vector< Number > has a bigger 
>> interface, PETScWrappers::MPI::Vector is extensively used in the examples. 
>> In which situations should we use each of them? Is there any known 
>> difference in performance? Thanks.
>>
>

-- 
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: PETScWrappers::MPI::Vector Vs parallel::distributed::Vector

[David F]

Hi, thanks for your answer. I have measured the time it takes for 
PETScWrappers::MPI::Vector, parallel::distributed::Vector< Number > and 
Vector< Number > to complete a very simple task, consisting in simple 
accessing the elements and assigning them to another variable, something 
like:

double a;
for(...)


On Friday, 26 August 2016 14:05:55 UTC+2, Bruno Turcksin wrote:
>
> Hi,
>
> I guess it's more a question of preference. What I do is using the same 
> vector type as the matrix type: PETSc matrix -> PETSc vector, Trilinos 
> matrix -> Trilinos vector, matrix-free -> deal.II vector. deal.II vector 
> can use multithreading unlike PETSc vector but if you are using MPI, I 
> don't think that you will see a big difference.
>
> Best,
>
> Bruno
>
> On Thursday, August 25, 2016 at 5:30:31 PM UTC-4, David F wrote:
>>
>> Hello, I would like to know if among PETScWrappers::MPI::Vector and 
>> parallel::distributed::Vector< Number >, one of them is preferred over the 
>> other. They both seem to have a similar functionality and a similar 
>> interface. Although parallel::distributed::Vector< Number > has a bigger 
>> interface, PETScWrappers::MPI::Vector is extensively used in the examples. 
>> In which situations should we use each of them? Is there any known 
>> difference in performance? Thanks.
>>
>

-- 
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] PETScWrappers::MPI::Vector Vs parallel::distributed::Vector

[David F]

Hello, I would like to know if among PETScWrappers::MPI::Vector and 
parallel::distributed::Vector< Number >, one of them is preferred over the 
other. They both seem to have a similar functionality and a similar 
interface. Although parallel::distributed::Vector< Number > has a bigger 
interface, PETScWrappers::MPI::Vector is extensively used in the examples. 
In which situations should we use each of them? Is there any known 
difference in performance? Thanks.

-- 
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.

Re: [deal.II] Efficient implementation of varying Dirichlet BCs when everything else in the system is constant

[David F]

Thank you, I'll take a look in detail at those examples. However, those 
examples don't use PETSc and the problem is that I can't find any method in 
PETScWrappers::SparseMatrix to copy a matrix. Which way would you suggest? 

David.


On Saturday, 20 August 2016 00:41:05 UTC+2, Wolfgang Bangerth wrote:
>
> On 08/19/2016 04:23 PM, David F wrote: 
> > 
> > P.S.: I tried to simply repeat the process to re-apply the BCs with a 
> new 
> > value using the same matrix, but I get always the first BCs I apply  (I 
> think 
> > because once I apply it, the matrix knows it has been already condensed 
> and 
> > ignores the next calls). 
>
> Yes. What you need to do is in every time step copy the matrix you have 
> previously assembled into a new matrix, call apply_boundary_values() to 
> this 
> matrix and your current right hand side, and then solve the linear system. 
> This way, you always start from the same matrix. 
>
> I think one of steps 23, 24, 25, 26 does this. 
>
> Best 
>   Wolfgang 
>
> -- 
>  
> 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 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] Efficient implementation of varying Dirichlet BCs when everything else in the system is constant

[David F]

Hello everyone, I looked in the forum, in the documentation and in the 
tutorials but I couldn't find a way to solve this problem in the way I 
need. I would appreciate a lot any help you could provide me. 


I would like to have Dirichlet BCs applied to faces with some boundary id. 
This id is not changing (i.e., the constrained faces are always the same), 
and the mesh, elastic properties etc. are also constant throughout the 
program. The function used to define these BCs has always the same shape, 
the only thing that is changing is the value of the Dirichlet BCs (e.g., 
imagine that first you constrain the displacement of a face with a 
spatially-constant value of 1.0 through the face, and in the next solution 
step to 2.0 and so on). I need to solve this system many times, each time 
with a different value of the Dirichlet BCs. To do this efficiently, I 
would like to precalculate as much as possible, for example, assembling the 
system matrix only once etc. 

Now my question is, given that the sparsity pattern is not changing (as the 
constrained nodes are the same), is there a way to solve this in a manner 
such that only the solution and rhs are changing from one step to the next? 
If the answer is that the system matrix must necessarily change, which 
should be the strategy to reassemble only the changing cells (and avoiding 
thus reassembling the whole matrix)? In summary, what's the most efficient 
way to repeatedly solve a system in which the only changes from step to 
step are the value (not the region over which they are applied) of the 
Dirichlet BCs?



About my current implementation:

I don't have hanging nodes, and I'm using PETScWrappers::MPI::
SparseMatrix and PETScWrappers::MPI::Vector.

I assemble the matrix with

void  distribute_local_to_global 

 
(const FullMatrix 
< typename 
MatrixType::value_type > &local_matrix, const std::vector< size_type 

 
> &local_dof_indices, MatrixType &global_matrix) const 
and apply the BCs with

void interpolate_boundary_values 

 
(const DoFHandlerType &dof, const types::boundary_id 

 
boundary_component, const Function 
< 
DoFHandlerType::space_dimension, double > &boundary_function, std::map< 
types::global_dof_index 
,
 
double > &boundary_values, const ComponentMask 
 
&component_mask=ComponentMask 
())

void apply_boundary_values 

 
(const std::map< types::global_dof_index 
,
 
double > &boundary_values, PETScWrappers::MPI::SparseMatrix 

 
&matrix, PETScWrappers::MPI::Vector 

 
&solution, PETScWrappers::MPI::Vector 

 
&right_hand_side, const bool eliminate_columns=true)




P.S.: I tried to simply repeat the process to re-apply the BCs with a new 
value using the same matrix, but I get always the first BCs I apply  (I 
think because once I apply it, the matrix knows it has been already 
condensed and ignores the next calls).

-- 
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] Efficient implementation of varying Dirichlet BCs when everything else in the system is constant

[David F]

Hello everyone, I looked in the forum, in the documentation and in the 
tutorials but I couldn't find a way to solve this problem in the way I 
need. I would appreciate a lot any help you could provide me. 


I would like to have Dirichlet BCs applied to faces with some boundary id. 
This id is not changing (i.e., the constrained faces are always the same), 
and the mesh, elastic properties etc. are also constant throughout the 
program. The function used to defined these BCs has always the same shape, 
the only thing that is changing is the value of the Dirichlet BCs (e.g., 
imagine that first you constrain the displacement of a face with a 
spatially-constant value of 1.0 through the face, and in the next solution 
step to 2.0 and so on). I need to solve this system many times, each time 
with a different value of the Dirichlet BCs. To do this efficiently, I 
would like to precalculate as much as possible, for example, assembling the 
system matrix only once etc. 

Now my question is, given that the sparsity pattern is not changing (as the 
constraints nodes are the same), is there a way to solve this in a manner 
such that only the solution and rhs are changing from one step to the next? 
If the answer is that the system matrix must necessarily change, which 
should be the strategy to reassemble only the changing cells (and avoiding 
thus reassembling the whole matrix)? In summary, what's the most efficient 
way to repeatedly solve a system in which the only changes from step to 
step are the value (not the region over which they are applied) of the 
Dirichlet BCs?



About my current implementation:

I don't have hanging nodes, and I'm using PETScWrappers::MPI::SparseMatrix 
and PETScWrappers::MPI::Vector.

I assemble the matrix with

void  distribute_local_to_global 

 
(const FullMatrix 
< typename 
MatrixType::value_type > &local_matrix, const std::vector< size_type 

 
> &local_dof_indices, MatrixType &global_matrix) const 
and apply the BCs with

void interpolate_boundary_values 

 
(const DoFHandlerType &dof, const types::boundary_id 

 
boundary_component, const Function 
< 
DoFHandlerType::space_dimension, double > &boundary_function, std::map< 
types::global_dof_index 
,
 
double > &boundary_values, const ComponentMask 
 
&component_mask=ComponentMask 
())

void apply_boundary_values 

 
(const std::map< types::global_dof_index 
,
 
double > &boundary_values, PETScWrappers::MPI::SparseMatrix 

 
&matrix, PETScWrappers::MPI::Vector 

 
&solution, PETScWrappers::MPI::Vector 

 
&right_hand_side, const bool eliminate_columns=true)




P.S.: I tried to simply repeat the process to re-apply the BCs with a new 
value using the same matrix, but I get always the first BCs I apply  (I 
think because once I apply it, the matrix knows it has been already 
condensed and ignores the next calls).

















-- 
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] Problem making Petsc solver wrappers members of a class

[David F]

Hi, I'm working with step-17 and I would like to extend it to fit to my 
simulation. The requirement is that the solver must be called over the same 
system matrix a big number of times. Only the RHS is changing, so I would 
like to have the solver_control, solver cg and preconditioner as members of 
the class, instead of creating them every iteration as in the example.

However, making any of these 3 objects a member of the class an 
constructing it in the initialization list of the class, gives different 
errors. I think this is a problem related with the Petsc solvers because I 
could do it with the normal iterative solvers without Petsc. Is there a way 
to avoid creating these 3 objects every time? Is it even worth it, or these 
objects are very light? Please take into account I must repeat it these 
step ~50,000 per simulation.

Thank you very much.
David.

-- 
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.