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

```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 means the resolution of the inhomogeneities. For n=1 each```

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

```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
> > 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
---
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 view this discussion on the web visit

```

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

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

David.

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
---
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 view this discussion on the web visit

```

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

```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
---
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 view this discussion on the web visit

```

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

```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
---
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 view this discussion on the web visit

```

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

```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
---
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 view this discussion on the web visit

```

### [deal.II] Periodicity with simultaneous displacement controlled BCs

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

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
---
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 view this discussion on the web visit

```

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

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

```

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

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

```

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

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

```

### [deal.II] Obtaining the final position of the vertices

```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?

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
---
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

```

### [deal.II] Calculating gradients at arbitrary points

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

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

```

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

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

```

### Re: [deal.II] matrix factorization reuse?

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

```

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

```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?

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 > _matrix, const std::vector< size_type

> _dof_indices, MatrixType _matrix) const
and apply the BCs with

void interpolate_boundary_values

(const DoFHandlerType , const types::boundary_id

boundary_component, const Function
<
DoFHandlerType::space_dimension, double > _function, std::map<
types::global_dof_index
,

())

void apply_boundary_values

(const std::map< types::global_dof_index
,

double > _values, PETScWrappers::MPI::SparseMatrix

, PETScWrappers::MPI::Vector

, PETScWrappers::MPI::Vector

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

```

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

```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?

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 > _matrix, const std::vector< size_type

> _dof_indices, MatrixType _matrix) const
and apply the BCs with

void interpolate_boundary_values

(const DoFHandlerType , const types::boundary_id

boundary_component, const Function
<
DoFHandlerType::space_dimension, double > _function, std::map<
types::global_dof_index
,

())

void apply_boundary_values

(const std::map< types::global_dof_index
,

double > _values, PETScWrappers::MPI::SparseMatrix

, PETScWrappers::MPI::Vector

, PETScWrappers::MPI::Vector

_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