[petsc-users] DMPlex with spring elements

Wed Sep 24 12:01:39 CDT 2014

On Wed, Sep 24, 2014 at 11:31 AM, Miguel Angel Salazar de Troya <
salazardetroya at gmail.com> wrote:

> Thanks for your response. My discretization is based on spring elements.
> For the linear one dimensional case in which each spring has a coefficient
> k, their jacobian would be this two by two matrix.
>
> [  k    -k ]
> [ -k     k ]
>
> and the internal force
>
> [ k ( Ui - Uj) ]
> [ k ( Uj - Ui) ]
>
> where Ui and Uj are the node displacements (just one displacement per node
> because it's one dimensional)
>
> For the two dimensional case, assuming small deformations, we have a
> four-by-four matrix. Each node has two degrees of freedom. We obtain it by
> performing the outer product of the vector (t , -t) where "t" is the vector
> that connects both nodes in a spring. This is for the case of small
> deformations. I would need to assemble each spring contribution to the
> jacobian and the residual like they were finite elements. The springs share
> nodes, that's how they are connected. This example is just the linear case,
> I will have to implement a nonlinear case in a similar fashion.
>
> Seeing the DMNetwork example, I think it's what I need, although I don't
> know much of power electric grids and it's hard for me to understand what's
> going on. Do you have a good reference to be able to follow the code? For
> example, why are they adding components to the edges?
>

Okay, I understand what you want now. I would recommend doing a simple
example by hand first. The PetscSection is very simple,
just 'dim' degrees of freedom for each vertex. It should look something like

DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
PetscSectionCreate(comm, &s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetChart(s, vStart, vEnd);
for (v = vStart; v < vEnd; ++v) {
  PetscSectionSetDof(s, v, dim);
}
PetscSectionSetUp(s);
DMSetDefaultSection(dm, s);
PetscSectionDestroy(&s);

Then when you form the residual, you loop over edges and add a contribution
to each vertex (I think, its still not
clear to me how your residual is defined)

VecGetArray(locF, &residual);
DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
for (c = cStart; c < cEnd; ++c) {
  const PetscInt *cone;
  PetscInt offA, offB;

  DMPlexGetCone(c, &cone);
  PetscSectionGetOffset(s, cone[0], &offA);
  PetscSectionGetOffset(s, cone[1], &offB);
  residual[offA] = <contribution for vertex A>;
  residual[offB] = <contribution for vertex B>;
}
VecRestoreArray(locF, &residual);
<do LocalToGlobal unless PETSc is doing it for you>

After that works, using DMNetwork should be much easier.

  Thanks,

    Matt

> 475:     DMNetworkAddComponent <http://www.mcs.anl.gov/petsc/petsc-as/petsc-current/docs/manualpages/DM/DMNetworkAddComponent.html#DMNetworkAddComponent>(networkdm,i,componentkey[0],&pfdata.branch[i-eStart]);
>
> Miguel
>
>
> On Tue, Sep 23, 2014 at 11:13 PM, Abhyankar, Shrirang G. <
> abhyshr at mcs.anl.gov> wrote:
>
>> You may also want to take a look at the DMNetwork framework that can be
>> used for general unstructured networks that don't use PDEs. Its
>> description is given in the manual and an example is in
>> src/snes/examples/tutorials/network/pflow.
>>
>> Shri
>>
>> From:  Matthew Knepley <knepley at gmail.com>
>> Date:  Tue, 23 Sep 2014 22:40:52 -0400
>> To:  Miguel Angel Salazar de Troya <salazardetroya at gmail.com>
>> Cc:  "petsc-users at mcs.anl.gov" <petsc-users at mcs.anl.gov>
>> Subject:  Re: [petsc-users] DMPlex with spring elements
>>
>>
>> >On Tue, Sep 23, 2014 at 4:01 PM, Miguel Angel Salazar de Troya
>> ><salazardetroya at gmail.com> wrote:
>> >
>> >Hi all
>> >I was wondering if it could be possible to build a model similar to the
>> >example snes/ex12.c, but with spring elements (for elasticity) instead of
>> >simplicial elements. Spring elements in a grid, therefore each element
>> >would have two nodes and each node two components. There would be more
>> >differences, because instead of calling the functions f0,f1,g0,g1,g2 and
>> >g3 to build the residual and the jacobian, I would call a routine that
>> >would build the residual vector and the jacobian matrix directly. I would
>> >not have shape functions whatsoever. My problem is discrete, I don't have
>> >a PDE and my equations are algebraic. What is the best way in petsc to
>> >solve this problem? Is there any example that I can follow? Thanks in
>> >advance
>> >
>> >
>> >
>> >
>> >Yes, ex12 is fairly specific to FEM. However, I think the right tools for
>> >what you want are
>> >DMPlex and PetscSection. Here is how I would proceed:
>> >
>> >  1) Make a DMPlex that encodes a simple network that you wish to
>> simulate
>> >
>> >  2) Make a PetscSection that gets the data layout right. Its hard from
>> >the above
>> >      for me to understand where you degrees of freedom actually are.
>> >This is usually
>> >      the hard part.
>> >
>> >  3) Calculate the residual, so you can check an exact solution. Here you
>> >use the
>> >      PetscSectionGetDof/Offset() for each mesh piece that you are
>> >interested in. Again,
>> >      its hard to be more specific when I do not understand your
>> >discretization.
>> >
>> >  Thanks,
>> >
>> >    Matt
>> >
>> >
>> >Miguel
>> >
>> >
>> >
>> >--
>> >Miguel Angel Salazar de Troya
>> >Graduate Research Assistant
>> >Department of Mechanical Science and Engineering
>> >University of Illinois at Urbana-Champaign
>> >(217) 550-2360 <tel:%28217%29%20550-2360>
>> >salaza11 at illinois.edu
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >--
>> >What most experimenters take for granted before they begin their
>> >experiments is infinitely more interesting than any results to which
>> >their experiments lead.
>> >-- Norbert Wiener
>>
>>
>
>
> --
> *Miguel Angel Salazar de Troya*
> Graduate Research Assistant
> Department of Mechanical Science and Engineering
> University of Illinois at Urbana-Champaign
> (217) 550-2360
> salaza11 at illinois.edu
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140924/eec95f25/attachment.html>