libmesh-users — General information for libmesh users

Re: [Libmesh-users] [RB] Question about a poseteriori error estimation

[<ss...@pu...>]

Thanks to your help, I understood inner-product and norm in the LibMesh.

 

However, I have another question about the dual norm of residual.

Despite the parameter-independent norm, I know that we need to a reference parameter set to compute the residual.

 

How is this reference parameter set selected in our parameter domain? 

Do they use random or middle value in the parameter domain as a reference parameter set?

 

I will be waiting for your reply.

 

Best regards,

SKang

 

 

From: David Knezevic <dav...@ak...> 
Sent: Monday, October 15, 2018 2:45 AM
To: 강신성 <ss...@pu...>
Cc: libmesh-users <lib...@li...>
Subject: Re: [Libmesh-users] [RB] Question about a poseteriori error estimation

 

On Sun, Oct 14, 2018 at 3:45 AM <ss...@pu... <mailto:ss...@pu...> > wrote:

Thank you for your prompt reply, David.

 

I looked at the “assembly.h” in RB example 1 and my code, and I understood your answer.

You said inner-product must be parameter-independent.

That is, dual norm of residual in LibMesh is always “X-norm” (parameter-independent norm), right?

 

Yes, it uses the parameter-independent inner product and norm that you specify. (I'm not familiar with referring to this as the "X-norm", but if that's what you want to call it then that's fine.)

 

By the way, if you do not mind, could you please tell me what special case used the parameter in “compute_residual_dual_norm(N)” code?

 

You call compute_residual_dual_norm after you do a solve. So the parameter mu that is used inside compute_residual_dual_norm is the parameter that was used for the preceding solve. The idea is that you're getting the dual norm of the residual associated with your previous solve.

 

Best,

David

 

 

 

From: David Knezevic <dav...@ak... <mailto:dav...@ak...> > 
Sent: Friday, October 12, 2018 9:30 PM
To: 강신성 <ss...@pu... <mailto:ss...@pu...> >
Cc: libmesh-users <lib...@li... <mailto:lib...@li...> >
Subject: Re: [Libmesh-users] [RB] Question about a poseteriori error estimation

 

On Fri, Oct 12, 2018 at 2:58 AM <ss...@pu... <mailto:ss...@pu...> > wrote:

Hello, all.



I have a question about the RB error estimation in LibMesh.



I know that the maximum error bound is computed in the function
"RBConstruction::compute_max_error_bound()."

Here we should compute the "dual norm of residual" using the function
"compute_residual_dual_norm(N)."

But I do not know exactly if this "dual norm of residual" is "X-norm" or
"energy-norm," where "X-norm" is parameter-independent and "energy-norm" is
parameter-dependent.

 

The norm that is used is specified by the user. For example, in reduced_basis_ex1, we set the inner product in SimpleRBConstruction to be the "A0" operator (and the inner product is used to generate the norm). So the inner-product and norm that is used is completely up to you.

 

Note, though, that the inner product must be parameter independent. This is necessary since we don't want to change the inner product when we change parameters.

 

 

When I looked at the "compute_residual_dual_norm(N)" code, the current
parameter is called by the "const RBParameters & mu = get_parameters();"
code line.

So I guess the "dual norm of residual" is the "energy-norm," and I would
like to ask if it is correct.

 

No, the residual is whatever is specified by the user by calling "set_inner_product_assembly". The parameter is used to evaluate the dual norm in any specific case, but that is not related to the inner product that is used.

 

Best,
David

Re: [Libmesh-users] [RB] Question about a poseteriori error estimation

[David K. <dav...@ak...>]

On Sun, Oct 14, 2018 at 3:45 AM <ss...@pu...> wrote:

> Thank you for your prompt reply, David.
>
>
>
> I looked at the “assembly.h” in RB example 1 and my code, and I
> understood your answer.
>
> You said inner-product must be parameter-independent.
>
> That is, dual norm of residual in LibMesh is always “X-norm” (parameter-independent
> norm), right?
>

Yes, it uses the parameter-independent inner product and norm that you
specify. (I'm not familiar with referring to this as the "X-norm", but if
that's what you want to call it then that's fine.)


> By the way, if you do not mind, could you please tell me what special
> case used the parameter in “compute_residual_dual_norm(N)” code?
>

You call compute_residual_dual_norm after you do a solve. So the parameter
mu that is used inside compute_residual_dual_norm is the parameter that was
used for the preceding solve. The idea is that you're getting the dual norm
of the residual associated with your previous solve.

Best,
David




> *From:* David Knezevic <dav...@ak...>
> *Sent:* Friday, October 12, 2018 9:30 PM
> *To:* 강신성 <ss...@pu...>
> *Cc:* libmesh-users <lib...@li...>
> *Subject:* Re: [Libmesh-users] [RB] Question about a poseteriori error
> estimation
>
>
>
> On Fri, Oct 12, 2018 at 2:58 AM <ss...@pu...> wrote:
>
> Hello, all.
>
>
>
> I have a question about the RB error estimation in LibMesh.
>
>
>
> I know that the maximum error bound is computed in the function
> "RBConstruction::compute_max_error_bound()."
>
> Here we should compute the "dual norm of residual" using the function
> "compute_residual_dual_norm(N)."
>
> But I do not know exactly if this "dual norm of residual" is "X-norm" or
> "energy-norm," where "X-norm" is parameter-independent and "energy-norm" is
> parameter-dependent.
>
>
>
> The norm that is used is specified by the user. For example, in
> reduced_basis_ex1, we set the inner product in SimpleRBConstruction to be
> the "A0" operator (and the inner product is used to generate the norm). So
> the inner-product and norm that is used is completely up to you.
>
>
>
> Note, though, that the inner product must be parameter independent. This
> is necessary since we don't want to change the inner product when we change
> parameters.
>
>
>
>
>
> When I looked at the "compute_residual_dual_norm(N)" code, the current
> parameter is called by the "const RBParameters & mu = get_parameters();"
> code line.
>
> So I guess the "dual norm of residual" is the "energy-norm," and I would
> like to ask if it is correct.
>
>
>
> No, the residual is whatever is specified by the user by calling
> "set_inner_product_assembly". The parameter is used to evaluate the dual
> norm in any specific case, but that is not related to the inner product
> that is used.
>
>
>
> Best,
> David
>

Re: [Libmesh-users] [RB] Question about a poseteriori error estimation

[<ss...@pu...>]

Thank you for your prompt reply, David.

 

I looked at the “assembly.h” in RB example 1 and my code, and I understood your answer.

You said inner-product must be parameter-independent.

That is, dual norm of residual in LibMesh is always “X-norm” (parameter-independent norm), right?

 

By the way, if you do not mind, could you please tell me what special case used the parameter in “compute_residual_dual_norm(N)” code?

 

I’ll be waiting for your reply.

 

Best regards,

SKang

 

From: David Knezevic <dav...@ak...> 
Sent: Friday, October 12, 2018 9:30 PM
To: 강신성 <ss...@pu...>
Cc: libmesh-users <lib...@li...>
Subject: Re: [Libmesh-users] [RB] Question about a poseteriori error estimation

 

On Fri, Oct 12, 2018 at 2:58 AM <ss...@pu... <mailto:ss...@pu...> > wrote:

Hello, all.



I have a question about the RB error estimation in LibMesh.



I know that the maximum error bound is computed in the function
"RBConstruction::compute_max_error_bound()."

Here we should compute the "dual norm of residual" using the function
"compute_residual_dual_norm(N)."

But I do not know exactly if this "dual norm of residual" is "X-norm" or
"energy-norm," where "X-norm" is parameter-independent and "energy-norm" is
parameter-dependent.

 

The norm that is used is specified by the user. For example, in reduced_basis_ex1, we set the inner product in SimpleRBConstruction to be the "A0" operator (and the inner product is used to generate the norm). So the inner-product and norm that is used is completely up to you.

 

Note, though, that the inner product must be parameter independent. This is necessary since we don't want to change the inner product when we change parameters.

 

 

When I looked at the "compute_residual_dual_norm(N)" code, the current
parameter is called by the "const RBParameters & mu = get_parameters();"
code line.

So I guess the "dual norm of residual" is the "energy-norm," and I would
like to ask if it is correct.

 

No, the residual is whatever is specified by the user by calling "set_inner_product_assembly". The parameter is used to evaluate the dual norm in any specific case, but that is not related to the inner product that is used.

 

Best,
David

Re: [Libmesh-users] [RB] Question about a poseteriori error estimation

[David K. <dav...@ak...>]

On Fri, Oct 12, 2018 at 2:58 AM <ss...@pu...> wrote:

> Hello, all.
>
>
>
> I have a question about the RB error estimation in LibMesh.
>
>
>
> I know that the maximum error bound is computed in the function
> "RBConstruction::compute_max_error_bound()."
>
> Here we should compute the "dual norm of residual" using the function
> "compute_residual_dual_norm(N)."
>
> But I do not know exactly if this "dual norm of residual" is "X-norm" or
> "energy-norm," where "X-norm" is parameter-independent and "energy-norm" is
> parameter-dependent.
>

The norm that is used is specified by the user. For example, in
reduced_basis_ex1, we set the inner product in SimpleRBConstruction to be
the "A0" operator (and the inner product is used to generate the norm). So
the inner-product and norm that is used is completely up to you.

Note, though, that the inner product must be parameter independent. This is
necessary since we don't want to change the inner product when we change
parameters.



> When I looked at the "compute_residual_dual_norm(N)" code, the current
> parameter is called by the "const RBParameters & mu = get_parameters();"
> code line.
>
> So I guess the "dual norm of residual" is the "energy-norm," and I would
> like to ask if it is correct.
>

No, the residual is whatever is specified by the user by calling
"set_inner_product_assembly". The parameter is used to evaluate the dual
norm in any specific case, but that is not related to the inner product
that is used.

Best,
David

[Libmesh-users] [RB] Question about a poseteriori error estimation

[<ss...@pu...>]

Hello, all.

 

I have a question about the RB error estimation in LibMesh.

 

I know that the maximum error bound is computed in the function
"RBConstruction::compute_max_error_bound()."

Here we should compute the "dual norm of residual" using the function
"compute_residual_dual_norm(N)."

But I do not know exactly if this "dual norm of residual" is "X-norm" or
"energy-norm," where "X-norm" is parameter-independent and "energy-norm" is
parameter-dependent.

 

When I looked at the "compute_residual_dual_norm(N)" code, the current
parameter is called by the "const RBParameters & mu = get_parameters();"
code line.

So I guess the "dual norm of residual" is the "energy-norm," and I would
like to ask if it is correct.

 

I always thank you for your help.

 

Regards,

SKang

Re: [Libmesh-users] Support for DGFEMContext

[David K. <dav...@ak...>]

On Thu, Oct 4, 2018 at 9:09 AM Stogner, Roy H <roy...@ic...>
wrote:

>
> On Wed, 3 Oct 2018, David Knezevic wrote:
>
> > On Wed, Oct 3, 2018 at 6:46 PM Salazar De Troya, Miguel
> > <sal...@ll...> wrote:
> >
> >> Would it be possible if you shared an example of how you use
> DGFEMContext?
> >
> > I don't have any example I can share easily, sorry.
> >
> > It's similar to how you use standard FEMContexts, but there is also space
> > for data from the neighbor element which is what you need in DG
> assembly...
>
> You've got a little bit of informative code in rb_construction.C, but
> none of the reduced_basis examples actually exercise the DG option.
> I'm sure you have a hundred higher priority items on your todo list,
> but even aside from the didactic benefits, it would be nice to get
> some test coverage there at some point.
>
> Paul has been working on a more arbitrary "FEMContext also initializes
> data from other elements" system for FSI codes, and when that's done
> I'll be tempted to refactor DGFEMContext to use it under the hood, and
> if that happens it'd be nice if you already had a "Roy can't break my
> stuff without noticing" safety net.
>

OK, good points.

I'll make myself a to-do to add a test case or example that uses
DGFEMContext.

Best,
David

Re: [Libmesh-users] Support for DGFEMContext

[Stogner, R. H <roy...@ic...>]

On Wed, 3 Oct 2018, David Knezevic wrote:

> On Wed, Oct 3, 2018 at 6:46 PM Salazar De Troya, Miguel
> <sal...@ll...> wrote:
>
>> Would it be possible if you shared an example of how you use DGFEMContext?
>
> I don't have any example I can share easily, sorry.
>
> It's similar to how you use standard FEMContexts, but there is also space
> for data from the neighbor element which is what you need in DG assembly...

You've got a little bit of informative code in rb_construction.C, but
none of the reduced_basis examples actually exercise the DG option.
I'm sure you have a hundred higher priority items on your todo list,
but even aside from the didactic benefits, it would be nice to get
some test coverage there at some point.

Paul has been working on a more arbitrary "FEMContext also initializes
data from other elements" system for FSI codes, and when that's done
I'll be tempted to refactor DGFEMContext to use it under the hood, and
if that happens it'd be nice if you already had a "Roy can't break my
stuff without noticing" safety net.
---
Roy

Re: [Libmesh-users] Support for DGFEMContext

[David K. <dav...@ak...>]

On Wed, Oct 3, 2018 at 6:46 PM Salazar De Troya, Miguel
<sal...@ll...> wrote:

Hi David,
>
>
>
> Would it be possible if you shared an example of how you use DGFEMContext?
>

I don't have any example I can share easily, sorry.

It's similar to how you use standard FEMContexts, but there is also space
for data from the neighbor element which is what you need in DG assembly...

Best,
David



>
>
> *From: *David Knezevic <dav...@ak...>
> *Date: *Monday, July 24, 2017 at 11:05 AM
> *To: *"Salazar De Troya, Miguel" <sal...@ll...>
> *Cc: *"lib...@li..." <
> lib...@li...>
> *Subject: *Re: [Libmesh-users] Support for DGFEMContext
>
>
>
> On Mon, Jul 24, 2017 at 1:24 PM, Salazar De Troya, Miguel <
> sal...@ll...> wrote:
>
> Hi
>
> I have seen this class could be useful for my problem, but I don’t see any
> examples using it. Does it work well? Has it been tested?
>
>
>
>
>
> Yes, it works, and I use it. I guess we never got around to adding an
> example or unit tests for it, though. If you feel like doing that, that
> would be appreciated :)
>
>
>
> David
>
>
>

Re: [Libmesh-users] Support for DGFEMContext

[Salazar De T. M. <sal...@ll...>]

Hi David,

Would it be possible if you shared an example of how you use DGFEMContext?

Thanks
Miguel

From: David Knezevic <dav...@ak...>
Date: Monday, July 24, 2017 at 11:05 AM
To: "Salazar De Troya, Miguel" <sal...@ll...>
Cc: "lib...@li..." <lib...@li...>
Subject: Re: [Libmesh-users] Support for DGFEMContext

On Mon, Jul 24, 2017 at 1:24 PM, Salazar De Troya, Miguel <sal...@ll...<mailto:sal...@ll...>> wrote:
Hi

I have seen this class could be useful for my problem, but I don’t see any examples using it. Does it work well? Has it been tested?


Yes, it works, and I use it. I guess we never got around to adding an example or unit tests for it, though. If you feel like doing that, that would be appreciated :)

David

[Libmesh-users] Job ad: FEA developer at Akselos in Boston

[David K. <dav...@ak...>]

Hi all,

At Akselos we're looking to hire an FEA developer in Boston, see the ad for
"FEA/Simulation Software Engineer" on this page
<https://akselos.com/careers/>.

Our platform uses libMesh extensively mostly for structural integrity
analysis applications, i.e. "Digital Twins". See here
<https://www.businesswire.com/news/home/20180927005503/en/Akselos-Innogy-Ventures-Shell-Ventures-Predictive-Digital>
or
our website <http://www.akselos.com> for more info about Akselos.

Best regards,
David

--

David J. Knezevic | CTO
Akselos | 210 Broadway, #201 | Cambridge, MA | 02139
Phone: +1-617-599-4755

 [image: www.akselos.com] <http://www.akselos.com>   [image:
@AkselosCAE] <http://www.twitter.com/akselosCAE>    [image:
linkedin.com/company/akselos]
<https://www.linkedin.com/company/akselos>

This e-mail and any attachments may contain confidential material for
the sole use of the intended recipient(s). Any review or distribution
by others is strictly prohibited. If you are not the intended
recipient, please contact the sender and delete all copies.

Re: [Libmesh-users] Quadrature rule for unit normals at centroids, edges, and nodes

[John P. <jwp...@gm...>]

On Mon, Sep 24, 2018 at 11:12 AM Amneet Bhalla <mai...@gm...>
wrote:

> Hi John,
>
> Is there a way to obtain phi values at the centroid of the element (or in
> general any other point inside the element)? We can potentially interpolate
> nodal normals to centroid normal.
>

QGauss at CONSTANT/FIRST order will give you a quadrature rule with one
point at the element centroid.

If you want values at an arbitrary point in the element, you can call the
FE::reinit() overload that takes pointers to vectors of points and weights
defining a "custom" quadrature rule.



>
> PPS: Just checking QTRAP is exact for FIRST order elements, and QSIMPSON is
> exact for SECOND order elements in libMesh.
>

Depends what you mean by "exact". QTrap/QSimpson will give you nodal
quadrature rules for FIRST and SECOND order elements, respectively. On a
per-node basis these rules are not as accurate as the corresponding Gauss
rules. For example, a 2x2 Gaussian quadrature rule can integrate a bi-cubic
function exactly, whereas the 2x2 QTrap rule should only be able to
integrate a bilinear function exactly IIRC.

-- 
John

Re: [Libmesh-users] Quadrature rule for unit normals at centroids, edges, and nodes

[Amneet B. <mai...@gm...>]

Hi John,

Is there a way to obtain phi values at the centroid of the element (or in
general any other point inside the element)? We can potentially interpolate
nodal normals to centroid normal.

For do not need exact values of the normal. It is just need to know if the
point is on one side of the element or the other for our application. If
obtaining phi at centroid is not easy, then we can simply average nodal
normals to get the centroid normal.

PPS: Just checking QTRAP is exact for FIRST order elements, and QSIMPSON is
exact for SECOND order elements in libMesh.

Thanks,
—Amneet

On Mon, Sep 24, 2018 at 7:38 AM John Peterson <jwp...@gm...> wrote:

> On Sat, Sep 22, 2018 at 8:42 AM Nishant Nangia <nis...@gm...
> >
> wrote:
>
> > Hi Folks,
> >
> > I am working on an application where I need the unit normals on all the
> > nodes, the centroid, and (in 3D) on the edges halfway between the nodes
> of
> > a surface mesh. The boilerplate code I am using to get the normals is
> > something like:
> >
> > libMesh::UniquePtr<FEBase> fe_bdry(FEBase::build(dim, fe_type));
> > libMesh::UniquePtr<QBase> qrule_bdry(fe_type.default_quadrature_rule(dim
> -
> > 1));
> > fe_bdry->attach_quadrature_rule(qrule_bdry.get());
> >
> > ...
> >
> > const std::vector<libMesh::Point>& bdry_normals = fe_bdry->get_normals();
> >
> > When constructing the QBase, is there a quadrature rule and order I can
> > pass to get the normals at those specific locations? My surface meshes
> will
> > always consist of EDGE2 elements in 2D and TRI3 elements in 3D.
> >
>
> The nodal quadrature rules in libmesh are
>
> QTrap (trapezoidal rule)
> QSimpson (Simpson's rule)
>
> but you should be aware that the normals you get at the nodes will depend
> on which element you are evaluating the normal from, that is, the element
> normal direction changes discontinuously at the nodes.
>
> --
> John
>
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>

Re: [Libmesh-users] Quadrature rule for unit normals at centroids, edges, and nodes

[John P. <jwp...@gm...>]

On Sat, Sep 22, 2018 at 8:42 AM Nishant Nangia <nis...@gm...>
wrote:

> Hi Folks,
>
> I am working on an application where I need the unit normals on all the
> nodes, the centroid, and (in 3D) on the edges halfway between the nodes of
> a surface mesh. The boilerplate code I am using to get the normals is
> something like:
>
> libMesh::UniquePtr<FEBase> fe_bdry(FEBase::build(dim, fe_type));
> libMesh::UniquePtr<QBase> qrule_bdry(fe_type.default_quadrature_rule(dim -
> 1));
> fe_bdry->attach_quadrature_rule(qrule_bdry.get());
>
> ...
>
> const std::vector<libMesh::Point>& bdry_normals = fe_bdry->get_normals();
>
> When constructing the QBase, is there a quadrature rule and order I can
> pass to get the normals at those specific locations? My surface meshes will
> always consist of EDGE2 elements in 2D and TRI3 elements in 3D.
>

The nodal quadrature rules in libmesh are

QTrap (trapezoidal rule)
QSimpson (Simpson's rule)

but you should be aware that the normals you get at the nodes will depend
on which element you are evaluating the normal from, that is, the element
normal direction changes discontinuously at the nodes.

-- 
John

[Libmesh-users] Quadrature rule for unit normals at centroids, edges, and nodes

[Nishant N. <nis...@gm...>]

Hi Folks,

I am working on an application where I need the unit normals on all the
nodes, the centroid, and (in 3D) on the edges halfway between the nodes of
a surface mesh. The boilerplate code I am using to get the normals is
something like:

libMesh::UniquePtr<FEBase> fe_bdry(FEBase::build(dim, fe_type));
libMesh::UniquePtr<QBase> qrule_bdry(fe_type.default_quadrature_rule(dim -
1));
fe_bdry->attach_quadrature_rule(qrule_bdry.get());

...

const std::vector<libMesh::Point>& bdry_normals = fe_bdry->get_normals();

When constructing the QBase, is there a quadrature rule and order I can
pass to get the normals at those specific locations? My surface meshes will
always consist of EDGE2 elements in 2D and TRI3 elements in 3D.

Thanks!
Nishant

*Nishant Nangia*
Northwestern University
Ph.D. Candidate | Engineering Sciences and Applied Mathematics
Tech L386

Re: [Libmesh-users] Build libmesh with VTK - VTK_INCLUDE?

[Yuxiang W. <yw...@vi...>]

Hi all,

Thank you so much for your help! John & Shayan - you are totally right that
I did not run `make install` after the `make`. Sorry about the stupid
mistake. I just started to learn more :)

And thank you Simone as well for the detailed explanation.

Thank you all again! It is exporting vtk perfectly now.

Best,
Shawn

On Thu, Sep 20, 2018 at 9:14 AM John Peterson <jwp...@gm...> wrote:

>
>
> On Wed, Sep 19, 2018 at 9:17 PM Yuxiang Wang <yw...@vi...> wrote:
>
>> Dear all,
>>
>> I'm trying to build libmesh with VTK and noticed that I need to provide
>> two
>> variables - the VTK_DIR and the VTK_INCLUDE.
>>
>> I downloaded the VTK source and built the library
>> with DVTK_Group_MPI:BOOL=ON, but could not find the include directory
>> after
>> I compiled VTK.
>
>
> After you compiled? Did you run "make install" yet?
>
> --
> John
>


-- 
Yuxiang "Shawn" Wang, PhD
yw...@vi...
+1 (434) 284-0836

Re: [Libmesh-users] metis and libmesh installed

[John P. <jwp...@gm...>]

On Thu, Sep 20, 2018 at 9:19 AM Mauro Cacace <ca...@gf...> wrote:

> Hi John,
>
> thanks for the quick reply.
>
> Please find attached the config.log.
>
> I will ask the IT admin for granting you access to the cluster.
>

Thanks.

In the meantime I can try and explain what we currently do to avoid this
issue. If you look at the value of the libmesh_optional_LIBS variable in
your config.log, you should see something like the following (newlines
added for readability):

-lglpk
-lz
-L/cluster/petsc/petsc-3.9.3/lib
-Wl,-rpath,/cluster/petsc/petsc-3.9.3/lib
-Wl,-rpath,/opt/intel/ipsxe_2017_update1/compilers_and_libraries_2017.1.132/linux/mpi/intel64/lib/debug_mt
-L/opt/intel/ipsxe_2017_update1/compilers_and_libraries_2017.1.132/linux/mpi/intel64/lib/debug_mt
-Wl,-rpath,/opt/intel/ipsxe_2017_update1/compilers_and_libraries_2017.1.132/linux/mpi/intel64/lib
-L/opt/intel/ipsxe_2017_update1/compilers_and_libraries_2017.1.132/linux/mpi/intel64/lib
-Wl,-rpath,/opt/gcc/gcc-7.2.0/lib/gcc/x86_64-pc-linux-gnu/7.2.0
-L/opt/gcc/gcc-7.2.0/lib/gcc/x86_64-pc-linux-gnu/7.2.0
-Wl,-rpath,/opt/gcc/gcc-7.2.0/lib64
-L/opt/gcc/gcc-7.2.0/lib64
-Wl,-rpath,/opt/intel/ipsxe_2017_update1/impi/2017.1.132/intel64/lib
-L/opt/intel/ipsxe_2017_update1/impi/2017.1.132/intel64/lib
-Wl,-rpath,/opt/intel/ipsxe_2017_update1/ipp/lib/intel64
-L/opt/intel/ipsxe_2017_update1/ipp/lib/intel64
-Wl,-rpath,/opt/intel/ipsxe_2017_update1/tbb/lib/intel64/gcc4.1
-L/opt/intel/ipsxe_2017_update1/tbb/lib/intel64/gcc4.1
-Wl,-rpath,/opt/intel/ipsxe_2017_update1/mkl/lib/intel64
-L/opt/intel/ipsxe_2017_update1/mkl/lib/intel64
-Wl,-rpath,/opt/intel/ipsxe_2017_update1/mkl/interfaces/fftw3xf
-L/opt/intel/ipsxe_2017_update1/mkl/interfaces/fftw3xf
-Wl,-rpath,/opt/gcc/gcc-7.2.0/lib
-L/opt/gcc/gcc-7.2.0/lib
-Wl,-rpath,/opt/intel/mpi-rt/2017.0.0/intel64/lib/debug_mt
-Wl,-rpath,/opt/intel/mpi-rt/2017.0.0/intel64/lib
-lpetsc
-lcmumps
-ldmumps
-lsmumps
-lzmumps
-lmumps_common
-lpord
-lscalapack
-lsuperlu_dist
-lHYPRE
-lspai
-lflapack
-lfblas
-lhdf5hl_fortran
-lhdf5_fortran
-lhdf5_hl
-lhdf5
-lparmetis
-lmetis
-lX11
-lmpifort
-lmpi
-lmpigi
-lrt
-lpthread
-lgfortran
-lm
-lgcc_s
-lquadmath
-lstdc++
-ldl
-Wl,-rpath,/usr/local/lib
-L/usr/local/lib


As you can see, "-L/usr/local/lib" appears *last* in this list, and we
actually do this on purpose to avoid linking in system libraries that may
be in conflict with custom-built software. I was under the impression that
you link against the *first* -lmetis that you come across, and that should
be the one in -L/cluster/petsc/petsc-3.9.3/lib, which comes first in your
link line. This impression may be incorrect, however.

One final comment: the error message you posted comes during the linking of
the meshid-opt executable, so it is possible that there is a libmesh_opt.so
in your libmesh/build/.libs directory? If so, I would be curious what the
output of "ldd" on that shared library is.

-- 
John

[Libmesh-users] Fwd: metis and libmesh installed

[John P. <jwp...@gm...>]

---------- Forwarded message ---------
From: John Peterson <jwp...@gm...>
Date: Thu, Sep 20, 2018 at 9:11 AM
Subject: Re: [Libmesh-users] metis and libmesh installed
To: <ca...@gf...>




On Thu, Sep 20, 2018 at 6:09 AM Mauro Cacace <ca...@gf...> wrote:

> Dear all,
>
> I am trying to install libmesh (as a part of the MOOSE framework) on a
> cluster (no root). The problem is that I am continuously having this
> error message:
>
> "/ld: gk_cur_jbufs: TLS definition
> in*/usr/local/lib/libmetis.a(error.c.o)* section .tdata mismatches
> non-TLS definition in /cluster/petsc/petsc-3.9.3/lib///libmetis.so
> <http://libmetis.so>//section .data//
> //*/cluster/petsc/petsc-3.9.3/lib/*/*/libmetis.so
> <http://libmetis.so>/*/: error adding symbols: Bad value//
> //make[1]: *** [meshid-opt] Error 1//
> //make[1]: *** Waiting for unfinished jobs..../"
>
> The problem is that there is a conflict of libraries (metis shared and
> static). The PETSc version installed (against which I would like to
> link) does provide a so library (under cluster/petsc/petsc-3.9.3/lib).
> However, libmesh always tries always to further link to an existing old
> static libraries in the common path (under /usr/local/lib/).
>
> I already asked on the MOOSE user group, and received this answer:
>
> /"//The issue seems to be that there are two libmetis libraries under
> consideration by the linker, and because their thread-local storage
> ("TLS") declarations disagree, the linker throws an error. The only way
> we have been able to fix this in the past is to remove the system
> libmetis (the one in /usr/local/lib) completely, but this is obviously
> not acceptable in all situations, so it would be nice to come up with a
> better solution. It would be nice, for example if the "first"
> //libmetis.so <http://libmetis.so>//could somehow take precedence and
> the second one could be ignored...//"/
>
> The IT administrator asks whether it is possible to exclude the static
> library from the making (I think it would not be possible to delete it).
> However, it is not clear to me if that is possible at all, and where.
>
> Any suggestion/help/way out would be really appreciated.
>

Hi Mauro,

Can you send me the config.log file from your build directory on this
machine (either direct email or Google drive link, this list does not allow
attachments).

Also, would it be possible for me to get temporary access to the cluster on
which you are experiencing the problem? None of our local machines exhibit
the issue, so it's impossible to reproduce/fix locally.

-- 
John


-- 
John

Re: [Libmesh-users] Build libmesh with VTK - VTK_INCLUDE?

[John P. <jwp...@gm...>]

On Wed, Sep 19, 2018 at 9:17 PM Yuxiang Wang <yw...@vi...> wrote:

> Dear all,
>
> I'm trying to build libmesh with VTK and noticed that I need to provide two
> variables - the VTK_DIR and the VTK_INCLUDE.
>
> I downloaded the VTK source and built the library
> with DVTK_Group_MPI:BOOL=ON, but could not find the include directory after
> I compiled VTK.


After you compiled? Did you run "make install" yet?

-- 
John

Re: [Libmesh-users] Build libmesh with VTK - VTK_INCLUDE?

[Simone R. <sim...@gm...>]

As Shayan said, you want to know where VTK is installed. If you have not specified the prefix when you installed VTK, you should be able to locate itineraries the terminal using “locate” (e.g. locate vtkCell.h). Note that this will also return the source VTK folder, which is not the installation folder.
Most likely the default prefix is /usr/local or /usr.

Once you know the install directory, say that the vtk prefix is VTKPREFIX  (e.g. VTKPREFIX=/usr/local/vtk), then, to configure libMesh with VTK 7.1.1, you can use the following flags
--enable-vtk  --with-vtk-include=$VTKPREFIX/include/vtk-7.1 --with-vtk-lib=$VTKPREFIX/lib

If you want to use a different VTK version you can change the number of the include folder.
I hope it helps,
Simone


> On Sep 20, 2018, at 1:28 AM, Shayan Hoshyari <s.h...@gm...> wrote:
> 
> When running cmake for VTK you must pass the install location (where the
> lib files and the header files will go) via
> -DCMAKE_INSTALL_PREFIX=/the/desired/path
> 
> Then running make install, will move everything to that location.
> 
> On Thu, Sep 20, 2018 at 5:17 AM Yuxiang Wang <yw...@vi...> wrote:
> 
>> Dear all,
>> 
>> I'm trying to build libmesh with VTK and noticed that I need to provide two
>> variables - the VTK_DIR and the VTK_INCLUDE.
>> 
>> I downloaded the VTK source and built the library
>> with DVTK_Group_MPI:BOOL=ON, but could not find the include directory after
>> I compiled VTK. Being curious, did anyone experienced the same thing
>> before? Would there be an other folder that I missed?
>> 
>> I am using Windows Subsystem Linux, Ubuntu 18.04.
>> 
>> Thank you,
>> Shawn
>> --
>> Yuxiang "Shawn" Wang, PhD
>> yw...@vi...
>> +1 (434) 284-0836
>> 
>> _______________________________________________
>> Libmesh-users mailing list
>> Lib...@li...
>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>> 
> -- 
> Shayan Hoshyari
> 
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users

[Libmesh-users] metis and libmesh installed

[Mauro C. <ca...@gf...>]

Dear all,

I am trying to install libmesh (as a part of the MOOSE framework) on a 
cluster (no root). The problem is that I am continuously having this 
error message:

"/ld: gk_cur_jbufs: TLS definition 
in*/usr/local/lib/libmetis.a(error.c.o)* section .tdata mismatches 
non-TLS definition in /cluster/petsc/petsc-3.9.3/lib///libmetis.so 
<http://libmetis.so>//section .data//
//*/cluster/petsc/petsc-3.9.3/lib/*/*/libmetis.so 
<http://libmetis.so>/*/: error adding symbols: Bad value//
//make[1]: *** [meshid-opt] Error 1//
//make[1]: *** Waiting for unfinished jobs..../"

The problem is that there is a conflict of libraries (metis shared and 
static). The PETSc version installed (against which I would like to 
link) does provide a so library (under cluster/petsc/petsc-3.9.3/lib). 
However, libmesh always tries always to further link to an existing old 
static libraries in the common path (under /usr/local/lib/).

I already asked on the MOOSE user group, and received this answer:

/"//The issue seems to be that there are two libmetis libraries under 
consideration by the linker, and because their thread-local storage 
("TLS") declarations disagree, the linker throws an error. The only way 
we have been able to fix this in the past is to remove the system 
libmetis (the one in /usr/local/lib) completely, but this is obviously 
not acceptable in all situations, so it would be nice to come up with a 
better solution. It would be nice, for example if the "first" 
//libmetis.so <http://libmetis.so>//could somehow take precedence and 
the second one could be ignored...//"/

The IT administrator asks whether it is possible to exclude the static 
library from the making (I think it would not be possible to delete it). 
However, it is not clear to me if that is possible at all, and where.

Any suggestion/help/way out would be really appreciated.

Thanks,

Mauro


-- 
Dr. Mauro Cacace
Department 6 Geotechnologies
Section 6.1 Basin Modelling
Phone : +49(0)331/288-1783
Fax : +49(0)331/288-1349
Email : ca...@gf...
_________________________________________

Helmholtz Centre Potsdam
*GFZ German Research Centre for Geosciences*
Foundation under public law of the federal state
of Brandenburg
Telegrafenberg, D-14473 Potsdam

Re: [Libmesh-users] Build libmesh with VTK - VTK_INCLUDE?

[Shayan H. <s.h...@gm...>]

When running cmake for VTK you must pass the install location (where the
lib files and the header files will go) via
-DCMAKE_INSTALL_PREFIX=/the/desired/path

Then running make install, will move everything to that location.

On Thu, Sep 20, 2018 at 5:17 AM Yuxiang Wang <yw...@vi...> wrote:

> Dear all,
>
> I'm trying to build libmesh with VTK and noticed that I need to provide two
> variables - the VTK_DIR and the VTK_INCLUDE.
>
> I downloaded the VTK source and built the library
> with DVTK_Group_MPI:BOOL=ON, but could not find the include directory after
> I compiled VTK. Being curious, did anyone experienced the same thing
> before? Would there be an other folder that I missed?
>
> I am using Windows Subsystem Linux, Ubuntu 18.04.
>
> Thank you,
> Shawn
> --
> Yuxiang "Shawn" Wang, PhD
> yw...@vi...
> +1 (434) 284-0836
>
> _______________________________________________
> Libmesh-users mailing list
> Lib...@li...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
-- 
Shayan Hoshyari

[Libmesh-users] Build libmesh with VTK - VTK_INCLUDE?

[Yuxiang W. <yw...@vi...>]

Dear all,

I'm trying to build libmesh with VTK and noticed that I need to provide two
variables - the VTK_DIR and the VTK_INCLUDE.

I downloaded the VTK source and built the library
with DVTK_Group_MPI:BOOL=ON, but could not find the include directory after
I compiled VTK. Being curious, did anyone experienced the same thing
before? Would there be an other folder that I missed?

I am using Windows Subsystem Linux, Ubuntu 18.04.

Thank you,
Shawn
-- 
Yuxiang "Shawn" Wang, PhD
yw...@vi...
+1 (434) 284-0836

Re: [Libmesh-users] Augmenting Sparsity to add Nonlocal DOFs

[Charlie T. <cha...@gm...>]

In particular, here is a slightly modified assemble function (from a
miscellaneous example for IPDG).  For some reason, this does not work

void assemble_nonlocal(EquationSystems & es,
                         const std::string & libmesh_dbg_var(system_name))
{
  libMesh::out << " assembling nonlocal stiffness... ";libMesh::out.flush();
  // // _NonlocalStiffness->zero();
  // EquationSystems & es = this->get_equation_systems();
  // const MeshBase & mesh = es.get_mesh();
  // const unsigned int dim = mesh.mesh_dimension();

  const MeshBase & mesh = es.get_mesh();
  const unsigned int dim = mesh.mesh_dimension();

  // Get a reference to the LinearImplicitSystem we are solving
  LinearImplicitSystem & nonlocal_system =
es.get_system<LinearImplicitSystem> ("NonlocalSystem");
  // Get some parameters that we need during assembly
  const Real penalty = es.parameters.get<Real> ("penalty");
  const unsigned int target_patch_size = es.parameters.get<unsigned int>
("target_patch_size");
  const Real kernelparam = es.parameters.get<Real> ("kernelparam");
  const Real horizon = es.parameters.get<Real> ("horizon");
  const Order _quad_order = es.parameters.get<Order> ("_quad_order");

  std::string refinement_type = es.parameters.get<std::string>
("refinement");

  // A reference to the DofMap object for this system.  The DofMap
  // object handles the index translation from node and element numbers
  // to degree of freedom numbers.  We will talk more about the DofMap
  const DofMap & dof_map = nonlocal_system.get_dof_map();

  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  FEType fe_type = nonlocal_system.variable_type(0);

  // Build a Finite Element object of the specified type.  Since the
  // FEBase::build() member dynamically creates memory we will
  // store the object as a std::unique_ptr<FEBase>.  This can be thought
  // of as a pointer that will clean up after itself.
  std::unique_ptr<FEBase> fe_elem  (FEBase::build(dim, fe_type));
  std::unique_ptr<FEBase> fe_neighbor  (FEBase::build(dim, fe_type));
  std::unique_ptr<FEBase> fe_elem_face(FEBase::build(dim, fe_type));
  std::unique_ptr<FEBase> fe_neighbor_face(FEBase::build(dim, fe_type));

  // Quadrature rules for numerical integration.
#ifdef _quad_type
  // QGauss qrule (dim, QORDER);
  // std::unique_ptr<QBase> qrule(QBase::build(_quad_type, dim,
_quad_order));
  std::unique_ptr<QBase> qrule(QBase::build(_quad_type, dim, _quad_order));
  std::unique_ptr<QBase> qrule_neigbor(QBase::build(_quad_type, dim,
_quad_order));
#else
  std::unique_ptr<QBase> qrule(QBase::build(0, dim,
fe_type.default_quadrature_order()));
  std::unique_ptr<QBase> qrule_neighbor(QBase::build(0, dim,
fe_type.default_quadrature_order()));
#endif
  // fe->attach_quadrature_rule (&qrule);
  fe_elem->attach_quadrature_rule (qrule.get());
  fe_neighbor->attach_quadrature_rule (qrule_neighbor.get());


#ifdef _quad_type
  std::unique_ptr<QBase> qface(QBase::build(_quad_type, dim, _quad_order));
#else
  std::unique_ptr<QBase> qface(QBase::build(0, dim,
fe_type.default_quadrature_order()));
#endif
  fe_elem_face->attach_quadrature_rule(qrule.get());
  fe_neighbor_face->attach_quadrature_rule(qface.get());







  const std::vector<std::vector<Real>> &  phi_elem = fe_elem->get_phi();
  const std::vector<std::vector<RealGradient>> & dphi_elem =
fe_elem->get_dphi();
  const std::vector<Real> & JxW_elem = fe_elem->get_JxW();
  const std::vector<Point> & qpoint_elem = fe_elem->get_xyz();

  const std::vector<std::vector<Real>> &  phi_neighbor =
fe_neighbor->get_phi();
  const std::vector<std::vector<RealGradient>> & dphi2 =
fe_neighbor->get_dphi();
  const std::vector<Real> & JxW_neighbor = fe_neighbor->get_JxW();
  const std::vector<Point> & qpoint_neighbor = fe_neighbor->get_xyz();

  const std::vector<Point> & qpoint_elem_global = fe_elem->get_xyz();
  const std::vector<Point> & qpoint_neighbor_global =
fe_neighbor->get_xyz();


  // Data for surface integrals on the element boundary
  const std::vector<std::vector<Real>> &  phi_face =
fe_elem_face->get_phi();
  const std::vector<std::vector<RealGradient>> & dphi_face =
fe_elem_face->get_dphi();
  const std::vector<Real> & JxW_face = fe_elem_face->get_JxW();
  const std::vector<Point> & qface_normals = fe_elem_face->get_normals();
  const std::vector<Point> & qface_points = fe_elem_face->get_xyz();

  // // Data for surface integrals on the neighbor boundary
  const std::vector<std::vector<Real>> &  phi_neighbor_face =
fe_neighbor_face->get_phi();
  const std::vector<std::vector<RealGradient>> & dphi_neighbor_face =
fe_neighbor_face->get_dphi();



  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;
  // Data structures to contain the element and neighbor boundary matrix
  // contribution. This matrices will do the coupling between the dofs of
  // the element and those of his neighbors.
  // Ken: matrix coupling elem and neighbor dofs
  DenseMatrix<Number> Kne;
  DenseMatrix<Number> Ken;
  DenseMatrix<Number> Kee;
  DenseMatrix<Number> Knn;


  DenseMatrix<Number> zero_matrix, zero_matrix1, zero_matrix2;

  std::vector<dof_id_type> dof_indices;
  libMesh::Patch::PMF patchtype = &Patch::add_face_neighbors;
  libMesh::Patch patch(mesh.processor_id());

    // for (const auto& elem : mesh.active_local_element_ptr_range())
    for (const auto& elem : mesh.active_element_ptr_range())
    {
      patch.build_around_element(elem, target_patch_size, patchtype);

      fe_elem->reinit(elem);
      dof_map.dof_indices (elem, dof_indices);
      const unsigned int n_dofs = dof_indices.size();
      Ke.resize (n_dofs, n_dofs);
      Fe.resize (n_dofs);
      Kee.resize (n_dofs, n_dofs);

      /* reposition */
      // Kee.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs,
n_u_dofs);
          libMesh::out << " iterate over elements in patch ";
libMesh::out.flush();

      for (auto neighbor : patch)
      {
        fe_neighbor->reinit(neighbor);
        std::vector<dof_id_type> neighbor_dof_indices;
        dof_map.dof_indices (neighbor, neighbor_dof_indices);
        const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();

        Ken.resize (n_dofs, n_neighbor_dofs);
        Kne.resize (n_neighbor_dofs, n_dofs);
        Ken.resize (n_neighbor_dofs, n_neighbor_dofs);

        for (unsigned int qp=0; qp<qrule->n_points(); qp++)
        {
          // this->get_matrix("NonlocalStiffness").add_matrix(zero_matrix,
dof_indices, dof_indices);
          // this->get_matrix("NonlocalStiffness").add_matrix(zero_matrix,
neighbor_dof_indices, neighbor_dof_indices);
          // this->get_matrix("NonlocalStiffness").add_matrix(zero_matrix,
neighbor_dof_indices, dof_indices);

          libMesh::out << " fill with zeros 1 "; libMesh::out.flush();
          zero_matrix1.resize (n_dofs, n_neighbor_dofs);
          nonlocal_system.matrix->add_matrix(zero_matrix1, dof_indices,
neighbor_dof_indices);
          //
nonlocal_system.matrix->get_matrix("NonlocalStiffness").close();

          libMesh::out << " fill with zeros 2 "; libMesh::out.flush();
          zero_matrix2.resize (n_neighbor_dofs,n_dofs);
          nonlocal_system.matrix->add_matrix(zero_matrix2,
neighbor_dof_indices, dof_indices);
          //
nonlocal_system.matrix->get_matrix("NonlocalStiffness").close();

          for (unsigned int qp2=0; qp2 < qrule->n_points(); qp2++)
            {
              Real kernelval = kernel(qpoint_elem_global[qp],
qpoint_neighbor_global[qp2], kernelparam);
              libMesh::out << " kernelval " << kernelval;
              libMesh::out.flush();
              for (unsigned int i=0; i<n_dofs; i++)
                for (unsigned int j=0; j<n_neighbor_dofs; j++)
                    Ken(i,j) += JxW_neighbor[qp2] * JxW_elem[qp] *
kernelval * (phi_neighbor[i][qp2] * phi_elem[j][qp]);

              for (unsigned int i=0; i<n_neighbor_dofs; i++)
                for (unsigned int j=0; j<n_dofs; j++)
                    Kne(i,j) += JxW_neighbor[qp2] * JxW_elem[qp] *
kernelval * (phi_neighbor[i][qp2] * phi_elem[j][qp]);

              for (unsigned int i=0; i<n_dofs; i++)
                for (unsigned int j=0; j<n_dofs; j++)
                    Kee(i,j) += JxW_neighbor[qp2] * JxW_elem[qp] *
kernelval * (phi_elem[i][qp2] * phi_elem[j][qp]);

              for (unsigned int i=0; i<n_neighbor_dofs; i++)
                for (unsigned int j=0; j<n_neighbor_dofs; j++)
                    Knn(i,j) += JxW_neighbor[qp2] * JxW_elem[qp] *
kernelval * (phi_neighbor[i][qp2] * phi_neighbor[j][qp]);
              }


            // We consider Dirichlet bc imposed via the interior penalty
method
            for (auto side : elem->side_index_range())
              {
                if (elem->neighbor_ptr(side) == nullptr)
                  {
                    // Pointer to the element face
                    fe_elem_face->reinit(elem, side);

                    std::unique_ptr<const Elem> elem_side
(elem->build_side_ptr(side));
                    // h element dimension to compute the interior penalty
penalty parameter
                    const unsigned int elem_b_order = static_cast<unsigned
int> (fe_elem_face->get_order());
                    const double h_elem =
elem->volume()/elem_side->volume() * 1./pow(elem_b_order, 2.);

                    for (unsigned int qp=0; qp<qface->n_points(); qp++)
                      {
                        Number bc_value = exact_solution(qface_points[qp],
es.parameters, "null", "void");
                        for (unsigned int i=0; i<n_dofs; i++)
                          {
                            // Matrix contribution
                            for (unsigned int j=0; j<n_dofs; j++)
                              {
                                // stability
                                Ke(i,j) += JxW_face[qp] * penalty/h_elem *
phi_face[i][qp] * phi_face[j][qp];

                                // consistency
                                Ke(i,j) -=
                                  JxW_face[qp] *
                                  (phi_face[i][qp] *
(dphi_face[j][qp]*qface_normals[qp]) +
                                   phi_face[j][qp] *
(dphi_face[i][qp]*qface_normals[qp]));
                              }

                            // RHS contributions
                            // stability
                            Fe(i) += JxW_face[qp] * bc_value *
penalty/h_elem * phi_face[i][qp];
                            // consistency
                            Fe(i) -= JxW_face[qp] * dphi_face[i][qp] *
(bc_value*qface_normals[qp]);
                          }
                      }
                  }
                }

            libMesh::out << " constrain nonlocal stiffness...
";libMesh::out.flush();

              dof_map.constrain_element_matrix(Kee, dof_indices,
dof_indices);
              dof_map.constrain_element_matrix(Ken, dof_indices,
neighbor_dof_indices);
              dof_map.constrain_element_matrix(Kne, neighbor_dof_indices,
dof_indices);
              dof_map.constrain_element_matrix(Knn, neighbor_dof_indices,
neighbor_dof_indices);

            libMesh::out << " add nonlocal stiffness...
";libMesh::out.flush();
              nonlocal_system.matrix->add_matrix(Kee, dof_indices,
dof_indices);
              nonlocal_system.matrix->add_matrix(Ken, dof_indices,
neighbor_dof_indices);
              nonlocal_system.matrix->add_matrix(Kne, neighbor_dof_indices,
dof_indices);
              nonlocal_system.matrix->add_matrix(Knn, neighbor_dof_indices,
neighbor_dof_indices);
              }//qp

            }//patch
          dof_map.constrain_element_vector(Fe, dof_indices);
          nonlocal_system.rhs->add_vector(Fe, dof_indices);
        // this->_NonlocalStiffness->add_matrix(Kee, dof_indices);
    }//elem

  libMesh::out << " assembled nonlocal stiffness... ";libMesh::out.flush();
}//fcn




On Tue, Sep 18, 2018 at 4:05 AM Charlie Talbot <cha...@gm...> wrote:

> After looking at example 9 again, I have to ask:
> Better yet: is there a way to augment the sparsity pattern without this
> function? I've read that I should
>
> "
>   To do this, you need to cast your matrix to a PetscMatrix, then do
>   petsc_matrix->mat() to get the PETSc Mat that you can pass to
> MatSetOption.
>
>   PetscMatrix<Number> * NonlocalStiffness_petsc =
> &(this->get_matrix("NonlocalStiffness"));
>   MatSetOption( NonlocalStiffness_petsc.mat() ,
> MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
> "
>
> Suppose I added a sparse matrix "TestMat" to the system.  How do I do
> this, explicitly?
> Thanks!
>
> On Tue, Sep 18, 2018 at 3:40 AM Charlie Talbot <cha...@gm...>
> wrote:
>
>> Hi, I modified the files in miscellaneous example 9 to augment the
>> sparsity pattern to include elements found using Patch, and I can't seem to
>> get this to work, what am I missing? Thanks!
>>
>> #ifndef AUGMENT_SPARSITY_NONLOCAL_H
>> #define AUGMENT_SPARSITY_NONLOCAL_H
>>
>> #include "libmesh/ghosting_functor.h"
>> #include "libmesh/mesh_base.h"
>> #include "libmesh/patch.h"
>> using libMesh::Elem;
>> using libMesh::GhostingFunctor;
>> using libMesh::MeshBase;
>> using libMesh::boundary_id_type;
>> using libMesh::processor_id_type;
>>
>> // Convenient typedef for a map for (element, side id) --> element
>> neighbor
>> typedef std::map<std::pair<const Elem *, unsigned char>, const Elem *>
>> ElementSideMap;
>>
>> // And the inverse map, but without sides in the key
>> typedef std::map<const Elem *, const Elem *> ElementMap;
>>
>> class AugmentSparsityNonlocal : public GhostingFunctor
>> {
>> private:
>>
>>   /**
>>    * The Mesh we're calculating on
>>    */
>>   MeshBase & _mesh;
>>
>>   /**
>>    * A map from (lower element ID, side ID) to matching upper element ID.
>> Here "lower"
>>    * and "upper" refer to the lower and upper (wrt +z direction) sides of
>> the crack in
>>    * our mesh.
>>    */
>>   ElementSideMap _lower_to_upper;
>>
>>   /**
>>    * The inverse (ignoring sides) of the above map.
>>    */
>>   ElementMap _upper_to_lower;
>>
>>   /**
>>    * Boundary IDs for the lower and upper faces of the "crack" in the
>> mesh.
>>    */
>>   boundary_id_type _crack_boundary_lower, _crack_boundary_upper;
>>
>>   /**
>>    * Make sure we've been initialized before use
>>    */
>>   bool _initialized;
>>
>> public:
>>   unsigned int _target_patch_size;
>>
>>   /**
>>    * Constructor.
>>    */
>>   AugmentSparsityNonlocal(MeshBase & mesh,
>>     unsigned int target_patch_size);
>>                              // boundary_id_type crack_boundary_lower,
>>                              // boundary_id_type crack_boundary_upper);
>>
>>   /**
>>    * @return a const reference to the lower-to-upper element ID map.
>>    */
>>   const ElementSideMap & get_lower_to_upper() const;
>>
>>   /**
>>    * User-defined function to augment the sparsity pattern.
>>    */
>>   virtual void operator() (const MeshBase::const_element_iterator &
>> range_begin,
>>                            const MeshBase::const_element_iterator &
>> range_end,
>>                            processor_id_type p,
>>                            map_type & coupled_elements) override;
>>
>>   /**
>>    * Rebuild the cached _lower_to_upper map whenever our Mesh has
>>    * changed.
>>    */
>>   virtual void mesh_reinit () override;
>>
>>   /**
>>    * Update the cached _lower_to_upper map whenever our Mesh has been
>>    * redistributed.  We'll be lazy and just recalculate from scratch.
>>    */
>>   virtual void redistribute () override
>>   { this->mesh_reinit(); }
>>
>>
>> };
>>
>> #endif
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> // local includes
>> #include "augment_sparsity_nonlocal.h"
>>
>> // libMesh includes
>> #include "libmesh/elem.h"
>>
>> using namespace libMesh;
>>
>> AugmentSparsityNonlocal::AugmentSparsityNonlocal(MeshBase & mesh,
>>   unsigned int target_patch_size):
>>                                                        //
>> boundary_id_type crack_boundary_lower,
>>                                                        //
>> boundary_id_type crack_boundary_upper) :
>>   _mesh(mesh),
>>   // _crack_boundary_lower(crack_boundary_lower),
>>   // _crack_boundary_upper(crack_boundary_upper),
>>   _initialized(false)
>> {
>>   this->mesh_reinit();
>>   this->_target_patch_size = target_patch_size;
>>
>> }
>>
>>
>> const ElementSideMap & AugmentSparsityNonlocal::get_lower_to_upper() const
>> {
>>   libmesh_assert(this->_initialized);
>>   return _lower_to_upper;
>> }
>>
>> void AugmentSparsityNonlocal::mesh_reinit ()
>> {
>>   this->_initialized = true;
>>
>>   // Loop over all elements (not just local elements) to make sure we find
>>   // "neighbor" elements on opposite sides of the crack.
>>
>>   // Map from (elem, side) to centroid
>>   std::map<std::pair<const Elem *, unsigned char>, Point> lower_centroids;
>>   std::map<std::pair<const Elem *, unsigned char>, Point> upper_centroids;
>>
>>   // for (const auto & elem : _mesh.active_element_ptr_range())
>>   //   for (auto side : elem->side_index_range())
>>   //     if (elem->neighbor_ptr(side) == nullptr)
>>   //       {
>>   //         if (_mesh.get_boundary_info().has_boundary_id(elem, side,
>> _crack_boundary_lower))
>>   //           {
>>   //             std::unique_ptr<const Elem> side_elem =
>> elem->build_side_ptr(side);
>>
>>   //             lower_centroids[std::make_pair(elem, side)] =
>> side_elem->centroid();
>>   //           }
>>
>>   //         if (_mesh.get_boundary_info().has_boundary_id(elem, side,
>> _crack_boundary_upper))
>>   //           {
>>   //             std::unique_ptr<const Elem> side_elem =
>> elem->build_side_ptr(side);
>>
>>   //             upper_centroids[std::make_pair(elem, side)] =
>> side_elem->centroid();
>>   //           }
>>   //       }
>>
>>   // If we're doing a reinit on a distributed mesh then we may not see
>>   // all the centroids, or even a matching number of centroids.
>>   // std::size_t n_lower_centroids = lower_centroids.size();
>>   // std::size_t n_upper_centroids = upper_centroids.size();
>>   // libmesh_assert(n_lower_centroids == n_upper_centroids);
>>
>>   // Clear _lower_to_upper. This map will be used for matrix assembly
>> later on.
>>   _lower_to_upper.clear();
>>
>>   // Clear _upper_to_lower. This map will be used for element ghosting
>>   // on distributed meshes, communication send_list construction in
>>   // parallel, and sparsity calculations
>>   _upper_to_lower.clear();
>>
>>   // We do an N^2 search to find elements with matching centroids. This
>> could be optimized,
>>   // e.g. by first sorting the centroids based on their (x,y,z) location.
>>   {
>>     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator it
>>    = lower_centroids.begin();
>>     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
>> it_end = lower_centroids.end();
>>     for ( ; it != it_end; ++it)
>>       {
>>         Point lower_centroid = it->second;
>>
>>         // find closest centroid in upper_centroids
>>         Real min_distance = std::numeric_limits<Real>::max();
>>
>>         std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
>> inner_it     = upper_centroids.begin();
>>         std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
>> inner_it_end = upper_centroids.end();
>>
>>         for ( ; inner_it != inner_it_end; ++inner_it)
>>           {
>>             Point upper_centroid = inner_it->second;
>>
>>             Real distance = (upper_centroid - lower_centroid).norm();
>>             if (distance < min_distance)
>>               {
>>                 min_distance = distance;
>>                 _lower_to_upper[it->first] = inner_it->first.first;
>>               }
>>           }
>>
>>         // For pairs with local elements, we should have found a
>>         // matching pair by now.
>>         const Elem * elem     = it->first.first;
>>         const Elem * neighbor = _lower_to_upper[it->first];
>>         if (min_distance < TOLERANCE)
>>           {
>>             // fill up the inverse map
>>             _upper_to_lower[neighbor] = elem;
>>           }
>>         else
>>           {
>>             libmesh_assert_not_equal_to(elem->processor_id(),
>> _mesh.processor_id());
>>             // This must have a false positive; a remote element would
>>             // have been closer.
>>             _lower_to_upper.erase(it->first);
>>           }
>>       }
>>
>>     // Let's make sure we didn't miss any upper elements either
>> // #ifndef NDEBUG
>> //     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
>> inner_it     = upper_centroids.begin();
>> //     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
>> inner_it_end = upper_centroids.end();
>>
>> //     for ( ; inner_it != inner_it_end; ++inner_it)
>> //       {
>> //         const Elem * neighbor = inner_it->first.first;
>> //         if (neighbor->processor_id() != _mesh.processor_id())
>> //           continue;
>> //         ElementMap::const_iterator utl_it =
>> //           _upper_to_lower.find(neighbor);
>> //         libmesh_assert(utl_it != _upper_to_lower.end());
>> //       }
>> // #endif
>>   }
>> }
>>
>> void AugmentSparsityNonlocal::operator()
>>   (const MeshBase::const_element_iterator & range_begin,
>>    const MeshBase::const_element_iterator & range_end,
>>    processor_id_type p,
>>    map_type & coupled_elements)
>> {
>>   libmesh_assert(this->_initialized);
>>
>>   const CouplingMatrix * const null_mat = nullptr;
>>
>>   for (const auto & elem : as_range(range_begin, range_end))
>>     {
>>       if (elem->processor_id() != p)
>>         coupled_elements.insert (std::make_pair(elem, null_mat));
>>
>>       for (auto side : elem->side_index_range())
>>         if (elem->neighbor_ptr(side) == nullptr)
>>           {
>>             ElementSideMap::const_iterator ltu_it =
>>               _lower_to_upper.find(std::make_pair(elem, side));
>>             if (ltu_it != _lower_to_upper.end())
>>               {
>>                 const Elem * neighbor = ltu_it->second;
>>                 if (neighbor->processor_id() != p)
>>                   coupled_elements.insert (std::make_pair(neighbor,
>> null_mat));
>>               }
>>           }
>>
>>       ElementMap::const_iterator utl_it =
>>         _upper_to_lower.find(elem);
>>       if (utl_it != _upper_to_lower.end())
>>         {
>>           const Elem * neighbor = utl_it->second;
>>           if (neighbor->processor_id() != p)
>>             coupled_elements.insert (std::make_pair(neighbor, null_mat));
>>         }
>>
>>     }
>>
>>
>> /*
>> Nonlocal augment
>> */
>>   // const CouplingMatrix * const null_mat = nullptr;
>>   libMesh::Patch::PMF patchtype = &Patch::add_face_neighbors;
>>   // #ifndef this->_target_patch_size
>>   // const unsigned int _target_patch_size=4;
>>   // #endif
>> /* build element patches */
>>   for (const auto & elem : as_range(range_begin, range_end))
>>     {
>>       libMesh::Patch patch(_mesh.processor_id());
>>       patch.build_around_element(elem, _target_patch_size, patchtype);
>>       std::vector<const Elem *> patchvec;
>>       patchvec.reserve(patch.size());
>>       Patch::const_iterator        patch_it  = patch.begin();
>>       const Patch::const_iterator  patch_end = patch.end();
>>       for (; patch_it != patch_end; ++patch_it)
>>       {
>>         const Elem * elem2 = *patch_it;
>>         coupled_elements.insert (std::make_pair(elem2, null_mat));
>>           // patchvec.push_back(elem2);
>>       }
>>     }
>> }
>>
>>
>>
>> // MatSetOption(this->get_matrix("NonlocalStiffness"),
>> MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
>> //     DoFRenumbering::Cuthill_McKee (dof_handler);
>> //     hanging_node_constraints.clear();
>> //
>>  DoFTools::make_hanging_node_constraints(dof_handler,hanging_node_constraints);
>>
>> //     hanging_node_constraints.close();
>>
>> //     DynamicSparsityPattern dsp(dof_handler.n_dofs(),
>> dof_handler.n_dofs());
>> //     DoFTools::make_sparsity_pattern(dof_handler, dsp);
>> //     hanging_node_constraints.condense(dsp);
>> //     sparsity_pattern.copy_from(dsp);
>>
>> //     stiffness_matrix.reinit(sparsity_pattern);
>> //     mass_matrix.reinit(sparsity_pattern);
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>

Re: [Libmesh-users] Augmenting Sparsity to add Nonlocal DOFs

[Charlie T. <cha...@gm...>]

After looking at example 9 again, I have to ask:
Better yet: is there a way to augment the sparsity pattern without this
function? I've read that I should

"
  To do this, you need to cast your matrix to a PetscMatrix, then do
  petsc_matrix->mat() to get the PETSc Mat that you can pass to
MatSetOption.

  PetscMatrix<Number> * NonlocalStiffness_petsc =
&(this->get_matrix("NonlocalStiffness"));
  MatSetOption( NonlocalStiffness_petsc.mat() ,
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
"

Suppose I added a sparse matrix "TestMat" to the system.  How do I do this,
explicitly?
Thanks!

On Tue, Sep 18, 2018 at 3:40 AM Charlie Talbot <cha...@gm...> wrote:

> Hi, I modified the files in miscellaneous example 9 to augment the
> sparsity pattern to include elements found using Patch, and I can't seem to
> get this to work, what am I missing? Thanks!
>
> #ifndef AUGMENT_SPARSITY_NONLOCAL_H
> #define AUGMENT_SPARSITY_NONLOCAL_H
>
> #include "libmesh/ghosting_functor.h"
> #include "libmesh/mesh_base.h"
> #include "libmesh/patch.h"
> using libMesh::Elem;
> using libMesh::GhostingFunctor;
> using libMesh::MeshBase;
> using libMesh::boundary_id_type;
> using libMesh::processor_id_type;
>
> // Convenient typedef for a map for (element, side id) --> element neighbor
> typedef std::map<std::pair<const Elem *, unsigned char>, const Elem *>
> ElementSideMap;
>
> // And the inverse map, but without sides in the key
> typedef std::map<const Elem *, const Elem *> ElementMap;
>
> class AugmentSparsityNonlocal : public GhostingFunctor
> {
> private:
>
>   /**
>    * The Mesh we're calculating on
>    */
>   MeshBase & _mesh;
>
>   /**
>    * A map from (lower element ID, side ID) to matching upper element ID.
> Here "lower"
>    * and "upper" refer to the lower and upper (wrt +z direction) sides of
> the crack in
>    * our mesh.
>    */
>   ElementSideMap _lower_to_upper;
>
>   /**
>    * The inverse (ignoring sides) of the above map.
>    */
>   ElementMap _upper_to_lower;
>
>   /**
>    * Boundary IDs for the lower and upper faces of the "crack" in the mesh.
>    */
>   boundary_id_type _crack_boundary_lower, _crack_boundary_upper;
>
>   /**
>    * Make sure we've been initialized before use
>    */
>   bool _initialized;
>
> public:
>   unsigned int _target_patch_size;
>
>   /**
>    * Constructor.
>    */
>   AugmentSparsityNonlocal(MeshBase & mesh,
>     unsigned int target_patch_size);
>                              // boundary_id_type crack_boundary_lower,
>                              // boundary_id_type crack_boundary_upper);
>
>   /**
>    * @return a const reference to the lower-to-upper element ID map.
>    */
>   const ElementSideMap & get_lower_to_upper() const;
>
>   /**
>    * User-defined function to augment the sparsity pattern.
>    */
>   virtual void operator() (const MeshBase::const_element_iterator &
> range_begin,
>                            const MeshBase::const_element_iterator &
> range_end,
>                            processor_id_type p,
>                            map_type & coupled_elements) override;
>
>   /**
>    * Rebuild the cached _lower_to_upper map whenever our Mesh has
>    * changed.
>    */
>   virtual void mesh_reinit () override;
>
>   /**
>    * Update the cached _lower_to_upper map whenever our Mesh has been
>    * redistributed.  We'll be lazy and just recalculate from scratch.
>    */
>   virtual void redistribute () override
>   { this->mesh_reinit(); }
>
>
> };
>
> #endif
>
>
>
>
>
>
>
>
>
>
>
>
> // local includes
> #include "augment_sparsity_nonlocal.h"
>
> // libMesh includes
> #include "libmesh/elem.h"
>
> using namespace libMesh;
>
> AugmentSparsityNonlocal::AugmentSparsityNonlocal(MeshBase & mesh,
>   unsigned int target_patch_size):
>                                                        // boundary_id_type
> crack_boundary_lower,
>                                                        // boundary_id_type
> crack_boundary_upper) :
>   _mesh(mesh),
>   // _crack_boundary_lower(crack_boundary_lower),
>   // _crack_boundary_upper(crack_boundary_upper),
>   _initialized(false)
> {
>   this->mesh_reinit();
>   this->_target_patch_size = target_patch_size;
>
> }
>
>
> const ElementSideMap & AugmentSparsityNonlocal::get_lower_to_upper() const
> {
>   libmesh_assert(this->_initialized);
>   return _lower_to_upper;
> }
>
> void AugmentSparsityNonlocal::mesh_reinit ()
> {
>   this->_initialized = true;
>
>   // Loop over all elements (not just local elements) to make sure we find
>   // "neighbor" elements on opposite sides of the crack.
>
>   // Map from (elem, side) to centroid
>   std::map<std::pair<const Elem *, unsigned char>, Point> lower_centroids;
>   std::map<std::pair<const Elem *, unsigned char>, Point> upper_centroids;
>
>   // for (const auto & elem : _mesh.active_element_ptr_range())
>   //   for (auto side : elem->side_index_range())
>   //     if (elem->neighbor_ptr(side) == nullptr)
>   //       {
>   //         if (_mesh.get_boundary_info().has_boundary_id(elem, side,
> _crack_boundary_lower))
>   //           {
>   //             std::unique_ptr<const Elem> side_elem =
> elem->build_side_ptr(side);
>
>   //             lower_centroids[std::make_pair(elem, side)] =
> side_elem->centroid();
>   //           }
>
>   //         if (_mesh.get_boundary_info().has_boundary_id(elem, side,
> _crack_boundary_upper))
>   //           {
>   //             std::unique_ptr<const Elem> side_elem =
> elem->build_side_ptr(side);
>
>   //             upper_centroids[std::make_pair(elem, side)] =
> side_elem->centroid();
>   //           }
>   //       }
>
>   // If we're doing a reinit on a distributed mesh then we may not see
>   // all the centroids, or even a matching number of centroids.
>   // std::size_t n_lower_centroids = lower_centroids.size();
>   // std::size_t n_upper_centroids = upper_centroids.size();
>   // libmesh_assert(n_lower_centroids == n_upper_centroids);
>
>   // Clear _lower_to_upper. This map will be used for matrix assembly
> later on.
>   _lower_to_upper.clear();
>
>   // Clear _upper_to_lower. This map will be used for element ghosting
>   // on distributed meshes, communication send_list construction in
>   // parallel, and sparsity calculations
>   _upper_to_lower.clear();
>
>   // We do an N^2 search to find elements with matching centroids. This
> could be optimized,
>   // e.g. by first sorting the centroids based on their (x,y,z) location.
>   {
>     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator it
>    = lower_centroids.begin();
>     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
> it_end = lower_centroids.end();
>     for ( ; it != it_end; ++it)
>       {
>         Point lower_centroid = it->second;
>
>         // find closest centroid in upper_centroids
>         Real min_distance = std::numeric_limits<Real>::max();
>
>         std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
> inner_it     = upper_centroids.begin();
>         std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
> inner_it_end = upper_centroids.end();
>
>         for ( ; inner_it != inner_it_end; ++inner_it)
>           {
>             Point upper_centroid = inner_it->second;
>
>             Real distance = (upper_centroid - lower_centroid).norm();
>             if (distance < min_distance)
>               {
>                 min_distance = distance;
>                 _lower_to_upper[it->first] = inner_it->first.first;
>               }
>           }
>
>         // For pairs with local elements, we should have found a
>         // matching pair by now.
>         const Elem * elem     = it->first.first;
>         const Elem * neighbor = _lower_to_upper[it->first];
>         if (min_distance < TOLERANCE)
>           {
>             // fill up the inverse map
>             _upper_to_lower[neighbor] = elem;
>           }
>         else
>           {
>             libmesh_assert_not_equal_to(elem->processor_id(),
> _mesh.processor_id());
>             // This must have a false positive; a remote element would
>             // have been closer.
>             _lower_to_upper.erase(it->first);
>           }
>       }
>
>     // Let's make sure we didn't miss any upper elements either
> // #ifndef NDEBUG
> //     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
> inner_it     = upper_centroids.begin();
> //     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
> inner_it_end = upper_centroids.end();
>
> //     for ( ; inner_it != inner_it_end; ++inner_it)
> //       {
> //         const Elem * neighbor = inner_it->first.first;
> //         if (neighbor->processor_id() != _mesh.processor_id())
> //           continue;
> //         ElementMap::const_iterator utl_it =
> //           _upper_to_lower.find(neighbor);
> //         libmesh_assert(utl_it != _upper_to_lower.end());
> //       }
> // #endif
>   }
> }
>
> void AugmentSparsityNonlocal::operator()
>   (const MeshBase::const_element_iterator & range_begin,
>    const MeshBase::const_element_iterator & range_end,
>    processor_id_type p,
>    map_type & coupled_elements)
> {
>   libmesh_assert(this->_initialized);
>
>   const CouplingMatrix * const null_mat = nullptr;
>
>   for (const auto & elem : as_range(range_begin, range_end))
>     {
>       if (elem->processor_id() != p)
>         coupled_elements.insert (std::make_pair(elem, null_mat));
>
>       for (auto side : elem->side_index_range())
>         if (elem->neighbor_ptr(side) == nullptr)
>           {
>             ElementSideMap::const_iterator ltu_it =
>               _lower_to_upper.find(std::make_pair(elem, side));
>             if (ltu_it != _lower_to_upper.end())
>               {
>                 const Elem * neighbor = ltu_it->second;
>                 if (neighbor->processor_id() != p)
>                   coupled_elements.insert (std::make_pair(neighbor,
> null_mat));
>               }
>           }
>
>       ElementMap::const_iterator utl_it =
>         _upper_to_lower.find(elem);
>       if (utl_it != _upper_to_lower.end())
>         {
>           const Elem * neighbor = utl_it->second;
>           if (neighbor->processor_id() != p)
>             coupled_elements.insert (std::make_pair(neighbor, null_mat));
>         }
>
>     }
>
>
> /*
> Nonlocal augment
> */
>   // const CouplingMatrix * const null_mat = nullptr;
>   libMesh::Patch::PMF patchtype = &Patch::add_face_neighbors;
>   // #ifndef this->_target_patch_size
>   // const unsigned int _target_patch_size=4;
>   // #endif
> /* build element patches */
>   for (const auto & elem : as_range(range_begin, range_end))
>     {
>       libMesh::Patch patch(_mesh.processor_id());
>       patch.build_around_element(elem, _target_patch_size, patchtype);
>       std::vector<const Elem *> patchvec;
>       patchvec.reserve(patch.size());
>       Patch::const_iterator        patch_it  = patch.begin();
>       const Patch::const_iterator  patch_end = patch.end();
>       for (; patch_it != patch_end; ++patch_it)
>       {
>         const Elem * elem2 = *patch_it;
>         coupled_elements.insert (std::make_pair(elem2, null_mat));
>           // patchvec.push_back(elem2);
>       }
>     }
> }
>
>
>
> // MatSetOption(this->get_matrix("NonlocalStiffness"),
> MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
> //     DoFRenumbering::Cuthill_McKee (dof_handler);
> //     hanging_node_constraints.clear();
> //
>  DoFTools::make_hanging_node_constraints(dof_handler,hanging_node_constraints);
>
> //     hanging_node_constraints.close();
>
> //     DynamicSparsityPattern dsp(dof_handler.n_dofs(),
> dof_handler.n_dofs());
> //     DoFTools::make_sparsity_pattern(dof_handler, dsp);
> //     hanging_node_constraints.condense(dsp);
> //     sparsity_pattern.copy_from(dsp);
>
> //     stiffness_matrix.reinit(sparsity_pattern);
> //     mass_matrix.reinit(sparsity_pattern);
>
>
>
>
>
>
>
>
>
>

[Libmesh-users] Augmenting Sparsity to add Nonlocal DOFs

[Charlie T. <cha...@gm...>]

Hi, I modified the files in miscellaneous example 9 to augment the sparsity
pattern to include elements found using Patch, and I can't seem to get this
to work, what am I missing? Thanks!

#ifndef AUGMENT_SPARSITY_NONLOCAL_H
#define AUGMENT_SPARSITY_NONLOCAL_H

#include "libmesh/ghosting_functor.h"
#include "libmesh/mesh_base.h"
#include "libmesh/patch.h"
using libMesh::Elem;
using libMesh::GhostingFunctor;
using libMesh::MeshBase;
using libMesh::boundary_id_type;
using libMesh::processor_id_type;

// Convenient typedef for a map for (element, side id) --> element neighbor
typedef std::map<std::pair<const Elem *, unsigned char>, const Elem *>
ElementSideMap;

// And the inverse map, but without sides in the key
typedef std::map<const Elem *, const Elem *> ElementMap;

class AugmentSparsityNonlocal : public GhostingFunctor
{
private:

  /**
   * The Mesh we're calculating on
   */
  MeshBase & _mesh;

  /**
   * A map from (lower element ID, side ID) to matching upper element ID.
Here "lower"
   * and "upper" refer to the lower and upper (wrt +z direction) sides of
the crack in
   * our mesh.
   */
  ElementSideMap _lower_to_upper;

  /**
   * The inverse (ignoring sides) of the above map.
   */
  ElementMap _upper_to_lower;

  /**
   * Boundary IDs for the lower and upper faces of the "crack" in the mesh.
   */
  boundary_id_type _crack_boundary_lower, _crack_boundary_upper;

  /**
   * Make sure we've been initialized before use
   */
  bool _initialized;

public:
  unsigned int _target_patch_size;

  /**
   * Constructor.
   */
  AugmentSparsityNonlocal(MeshBase & mesh,
    unsigned int target_patch_size);
                             // boundary_id_type crack_boundary_lower,
                             // boundary_id_type crack_boundary_upper);

  /**
   * @return a const reference to the lower-to-upper element ID map.
   */
  const ElementSideMap & get_lower_to_upper() const;

  /**
   * User-defined function to augment the sparsity pattern.
   */
  virtual void operator() (const MeshBase::const_element_iterator &
range_begin,
                           const MeshBase::const_element_iterator &
range_end,
                           processor_id_type p,
                           map_type & coupled_elements) override;

  /**
   * Rebuild the cached _lower_to_upper map whenever our Mesh has
   * changed.
   */
  virtual void mesh_reinit () override;

  /**
   * Update the cached _lower_to_upper map whenever our Mesh has been
   * redistributed.  We'll be lazy and just recalculate from scratch.
   */
  virtual void redistribute () override
  { this->mesh_reinit(); }


};

#endif












// local includes
#include "augment_sparsity_nonlocal.h"

// libMesh includes
#include "libmesh/elem.h"

using namespace libMesh;

AugmentSparsityNonlocal::AugmentSparsityNonlocal(MeshBase & mesh,
  unsigned int target_patch_size):
                                                       // boundary_id_type
crack_boundary_lower,
                                                       // boundary_id_type
crack_boundary_upper) :
  _mesh(mesh),
  // _crack_boundary_lower(crack_boundary_lower),
  // _crack_boundary_upper(crack_boundary_upper),
  _initialized(false)
{
  this->mesh_reinit();
  this->_target_patch_size = target_patch_size;

}


const ElementSideMap & AugmentSparsityNonlocal::get_lower_to_upper() const
{
  libmesh_assert(this->_initialized);
  return _lower_to_upper;
}

void AugmentSparsityNonlocal::mesh_reinit ()
{
  this->_initialized = true;

  // Loop over all elements (not just local elements) to make sure we find
  // "neighbor" elements on opposite sides of the crack.

  // Map from (elem, side) to centroid
  std::map<std::pair<const Elem *, unsigned char>, Point> lower_centroids;
  std::map<std::pair<const Elem *, unsigned char>, Point> upper_centroids;

  // for (const auto & elem : _mesh.active_element_ptr_range())
  //   for (auto side : elem->side_index_range())
  //     if (elem->neighbor_ptr(side) == nullptr)
  //       {
  //         if (_mesh.get_boundary_info().has_boundary_id(elem, side,
_crack_boundary_lower))
  //           {
  //             std::unique_ptr<const Elem> side_elem =
elem->build_side_ptr(side);

  //             lower_centroids[std::make_pair(elem, side)] =
side_elem->centroid();
  //           }

  //         if (_mesh.get_boundary_info().has_boundary_id(elem, side,
_crack_boundary_upper))
  //           {
  //             std::unique_ptr<const Elem> side_elem =
elem->build_side_ptr(side);

  //             upper_centroids[std::make_pair(elem, side)] =
side_elem->centroid();
  //           }
  //       }

  // If we're doing a reinit on a distributed mesh then we may not see
  // all the centroids, or even a matching number of centroids.
  // std::size_t n_lower_centroids = lower_centroids.size();
  // std::size_t n_upper_centroids = upper_centroids.size();
  // libmesh_assert(n_lower_centroids == n_upper_centroids);

  // Clear _lower_to_upper. This map will be used for matrix assembly later
on.
  _lower_to_upper.clear();

  // Clear _upper_to_lower. This map will be used for element ghosting
  // on distributed meshes, communication send_list construction in
  // parallel, and sparsity calculations
  _upper_to_lower.clear();

  // We do an N^2 search to find elements with matching centroids. This
could be optimized,
  // e.g. by first sorting the centroids based on their (x,y,z) location.
  {
    std::map<std::pair<const Elem *, unsigned char>, Point>::iterator it
 = lower_centroids.begin();
    std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
it_end = lower_centroids.end();
    for ( ; it != it_end; ++it)
      {
        Point lower_centroid = it->second;

        // find closest centroid in upper_centroids
        Real min_distance = std::numeric_limits<Real>::max();

        std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
inner_it     = upper_centroids.begin();
        std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
inner_it_end = upper_centroids.end();

        for ( ; inner_it != inner_it_end; ++inner_it)
          {
            Point upper_centroid = inner_it->second;

            Real distance = (upper_centroid - lower_centroid).norm();
            if (distance < min_distance)
              {
                min_distance = distance;
                _lower_to_upper[it->first] = inner_it->first.first;
              }
          }

        // For pairs with local elements, we should have found a
        // matching pair by now.
        const Elem * elem     = it->first.first;
        const Elem * neighbor = _lower_to_upper[it->first];
        if (min_distance < TOLERANCE)
          {
            // fill up the inverse map
            _upper_to_lower[neighbor] = elem;
          }
        else
          {
            libmesh_assert_not_equal_to(elem->processor_id(),
_mesh.processor_id());
            // This must have a false positive; a remote element would
            // have been closer.
            _lower_to_upper.erase(it->first);
          }
      }

    // Let's make sure we didn't miss any upper elements either
// #ifndef NDEBUG
//     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
inner_it     = upper_centroids.begin();
//     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
inner_it_end = upper_centroids.end();

//     for ( ; inner_it != inner_it_end; ++inner_it)
//       {
//         const Elem * neighbor = inner_it->first.first;
//         if (neighbor->processor_id() != _mesh.processor_id())
//           continue;
//         ElementMap::const_iterator utl_it =
//           _upper_to_lower.find(neighbor);
//         libmesh_assert(utl_it != _upper_to_lower.end());
//       }
// #endif
  }
}

void AugmentSparsityNonlocal::operator()
  (const MeshBase::const_element_iterator & range_begin,
   const MeshBase::const_element_iterator & range_end,
   processor_id_type p,
   map_type & coupled_elements)
{
  libmesh_assert(this->_initialized);

  const CouplingMatrix * const null_mat = nullptr;

  for (const auto & elem : as_range(range_begin, range_end))
    {
      if (elem->processor_id() != p)
        coupled_elements.insert (std::make_pair(elem, null_mat));

      for (auto side : elem->side_index_range())
        if (elem->neighbor_ptr(side) == nullptr)
          {
            ElementSideMap::const_iterator ltu_it =
              _lower_to_upper.find(std::make_pair(elem, side));
            if (ltu_it != _lower_to_upper.end())
              {
                const Elem * neighbor = ltu_it->second;
                if (neighbor->processor_id() != p)
                  coupled_elements.insert (std::make_pair(neighbor,
null_mat));
              }
          }

      ElementMap::const_iterator utl_it =
        _upper_to_lower.find(elem);
      if (utl_it != _upper_to_lower.end())
        {
          const Elem * neighbor = utl_it->second;
          if (neighbor->processor_id() != p)
            coupled_elements.insert (std::make_pair(neighbor, null_mat));
        }

    }


/*
Nonlocal augment
*/
  // const CouplingMatrix * const null_mat = nullptr;
  libMesh::Patch::PMF patchtype = &Patch::add_face_neighbors;
  // #ifndef this->_target_patch_size
  // const unsigned int _target_patch_size=4;
  // #endif
/* build element patches */
  for (const auto & elem : as_range(range_begin, range_end))
    {
      libMesh::Patch patch(_mesh.processor_id());
      patch.build_around_element(elem, _target_patch_size, patchtype);
      std::vector<const Elem *> patchvec;
      patchvec.reserve(patch.size());
      Patch::const_iterator        patch_it  = patch.begin();
      const Patch::const_iterator  patch_end = patch.end();
      for (; patch_it != patch_end; ++patch_it)
      {
        const Elem * elem2 = *patch_it;
        coupled_elements.insert (std::make_pair(elem2, null_mat));
          // patchvec.push_back(elem2);
      }
    }
}



// MatSetOption(this->get_matrix("NonlocalStiffness"),
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
//     DoFRenumbering::Cuthill_McKee (dof_handler);
//     hanging_node_constraints.clear();
//
 DoFTools::make_hanging_node_constraints(dof_handler,hanging_node_constraints);

//     hanging_node_constraints.close();

//     DynamicSparsityPattern dsp(dof_handler.n_dofs(),
dof_handler.n_dofs());
//     DoFTools::make_sparsity_pattern(dof_handler, dsp);
//     hanging_node_constraints.condense(dsp);
//     sparsity_pattern.copy_from(dsp);

//     stiffness_matrix.reinit(sparsity_pattern);
//     mass_matrix.reinit(sparsity_pattern);