Re: [deal.II] Write to std::cout using MPI from subroutine
Thanks! I changed it to dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0 and now it's working. I appreciate the help Thanks, Alex On Wednesday, April 24, 2024 at 3:32:36 PM UTC-4 Wolfgang Bangerth wrote: > > > if (dealii::Utilities::MPI::this_mpi_process == 0) > > [...] > > warning: the address of ‘unsigned int > > dealii::Utilities::MPI::this_mpi_process(MPI_Comm)’ will never be NULL > > You need to say > Utilities::MPI::this_mpi_process() > instead of > Utilities::MPI::this_mpi_process > The former calls the function. The latter takes the address of the > function. > > Best > W. > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/79f0dd55-841b-46f5-a48c-f99ab2ebae9bn%40googlegroups.com.
[deal.II] Write to std::cout using MPI from subroutine
Hi everyone, I am looking to replicate the functionality of a ConditionalOstream in one of my subroutines. In my main file, I use the stream "pcout" based on Step-17. However, pcout is not passed in to my subroutine. I have tried to guard my code using something like: if (dealii::Utilities::MPI::this_mpi_process == 0) std::cout << "Text" << std::endl; This is the approach that is used in the ConditionalOstream pcout. However, the text never actually prints, and instead I'm given a warning when I compile: warning: the address of ‘unsigned int dealii::Utilities::MPI::this_mpi_process(MPI_Comm)’ will never be NULL [-Waddress] I'm unclear why the mpi_process will never be NULL; it works fine in the rest of my code, just not in an external file. Can you think of an effective way to print for only one process that won't cause this warning? Best regards, Alex -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/553dbdd9-d8ca-4ea4-8393-79a11302e053n%40googlegroups.com.
Re: [deal.II] Importing nodal BCs and accessing vertices
Thanks, Luca. I will try to implement your suggestions.
Best regards,
Alex
On Thursday, April 11, 2024 at 4:03:08 PM UTC-4 luca@gmail.com wrote:
> If you have access to the DG, you can call
> dh.get_triangulation().n_active_cells().
>
> Take a look at the tests subdirectory to find more inspiration on more
> ways to exploit your information.
>
> For example, I’d use DoFTools::map_dofs_to_support_points, use that to
> generate an r tree of indices, and then just ask for the one nearest
> neighbor (in the tests subdirectory there are a few examples that show how
> to do this for collections of points).
>
> Luca
>
> Il giorno 11 apr 2024, alle ore 15:43, Alex Quinlan
> ha scritto:
>
>
>
> Hi Wolfgang,
>
> I've gotten an R-tree subroutine set up for applying my boundary
> conditions based on rtree_8h.html
> <https://www.dealii.org/current/doxygen/deal.II/rtree_8h.html>. It seems
> to work, but I have a couple difficulties that I've highlighted in red.
> Hopefully, others might find use in my adaptation of the example and in my
> questions.
>
>
> My subroutine takes the following arguments: the DOF handler, the
> constraints object, and a vector of a custom nBC structure for holding the
> nodal boundary conditions that I've read from the file.
>
> template
> void nodal_bcs_rtree(std::vector & nBCs,
> dealii::DoFHandler &
> dof_handler,
> dealii::AffineConstraints & constraints )
> {
>namespace bgi = boost::geometry::index;
>
>// Create the R-tree ---
>std::vector> all_boxes(*10*); // TODO Need a way
> to set the correct size for this
>std::vector false> >> cell_list(*10*);
>
>unsigned int i = 0;
>for (const auto : dof_handler.active_cell_iterators())
> if (cell->is_locally_owned())
> {
>all_boxes[i] = cell->bounding_box(); // prepare the bounding boxes to
> be packed into the r-tree
>cell_list[i++] = cell;// create a corresponding vector
> of cells accessors
> }
> const auto idx_tree = pack_rtree_of_indices(all_boxes);
>
>// Search the R-tree
>
>const double nrad = 1e-5; // Tolerance on vertex position TODO make
> relative to cell size
>
>for (nBC bc : nBCs)
> {
> std::vector closest_cell_idx;
> idx_tree.query(bgi::nearest(bc.coords,1),
> std::back_inserter(closest_cell_idx)); // bc.coords is type Point<3>
>
> if (closest_cell_idx[0] < i-1) // Check that the cell index number
> will correspond to an actual cell in cell_list
> {
> const auto cell = cell_list[closest_cell_idx[0]];
>
> for (unsigned int v = 0; vn_vertices() ; v++) {
>if (cell->vertex(v).distance(bc.coords) < nrad){
> ~~~Apply constraints ~~~
> }
> }
> }
> }
> }
>
>
>
> 1) I'm having trouble getting the number of cells, and then sizing my
> vectors appropriately. From what I can see, the DoFHandler won't provide
> the number of cells, only the number of DoFs. In the short term, I've just
> thrown in the number 100,000 as the vector size. In the example (
> https://www.dealii.org/current/doxygen/deal.II/rtree_8h.html ), the line:
>
>std::vector> all_boxes(tria
> <https://www.dealii.org/current/doxygen/deal.II/tria__description_8cc.html#ab695339a227bd369fcea767631247cff>
> .n_locally_owned_active_cells());
>
> gets used to set the vector size. From within my subroutine, I don't have
> access to the triangulation object. Am I missing a way to get
> "n_locally_owned_active_cells" from the DoFHandler? Or should I just pass
> this value in as an argument to my subroutine?
>
>
> 2) Does my approach of creating vectors to hold the bounding boxes and the
> cell pointers make sense? Is there a better way to get access to the cell
> once I've determined the BoundingBox index from the R-tree?
>
> 3) Somewhat related to question 1, is there a way to check that the
> 'closest_cell' was actually found? I had an issue with MPI where the
> target point 'bc.coords' was outside the bounding box of some processes,
> resulting in some errors. I set this (closest_cell_idx[0] < i-1) to
> check that an actual BoundingBox was found.
>
>
> Thanks for your suggestion on the R-tree. I'd welcome any critiques of my
> implementation.
>
> Cheers,
> Alex
> On Wednesday, November 15, 2023 at 8:35:49 AM UTC-5 Alex Quinlan wrote:
>
>> Thanks, Wolfgang. I will abide by those guidelines.
>>
Re: [deal.II] Importing nodal BCs and accessing vertices
Hi Wolfgang,
I've gotten an R-tree subroutine set up for applying my boundary conditions
based on rtree_8h.html
<https://www.dealii.org/current/doxygen/deal.II/rtree_8h.html>. It seems
to work, but I have a couple difficulties that I've highlighted in red.
Hopefully, others might find use in my adaptation of the example and in my
questions.
My subroutine takes the following arguments: the DOF handler, the
constraints object, and a vector of a custom nBC structure for holding the
nodal boundary conditions that I've read from the file.
template
void nodal_bcs_rtree(std::vector & nBCs,
dealii::DoFHandler & dof_handler,
dealii::AffineConstraints & constraints )
{
namespace bgi = boost::geometry::index;
// Create the R-tree ---
std::vector> all_boxes(*10*); // TODO Need a way to
set the correct size for this
std::vector >> cell_list(*10*);
unsigned int i = 0;
for (const auto : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
all_boxes[i] = cell->bounding_box(); // prepare the bounding boxes to
be packed into the r-tree
cell_list[i++] = cell;// create a corresponding vector
of cells accessors
}
const auto idx_tree = pack_rtree_of_indices(all_boxes);
// Search the R-tree
const double nrad = 1e-5; // Tolerance on vertex position TODO make
relative to cell size
for (nBC bc : nBCs)
{
std::vector closest_cell_idx;
idx_tree.query(bgi::nearest(bc.coords,1),
std::back_inserter(closest_cell_idx)); // bc.coords is type Point<3>
if (closest_cell_idx[0] < i-1) // Check that the cell index number
will correspond to an actual cell in cell_list
{
const auto cell = cell_list[closest_cell_idx[0]];
for (unsigned int v = 0; vn_vertices() ; v++) {
if (cell->vertex(v).distance(bc.coords) < nrad){
~~~Apply constraints ~~~
}
}
}
}
}
1) I'm having trouble getting the number of cells, and then sizing my
vectors appropriately. From what I can see, the DoFHandler won't provide
the number of cells, only the number of DoFs. In the short term, I've just
thrown in the number 100,000 as the vector size. In the example (
https://www.dealii.org/current/doxygen/deal.II/rtree_8h.html
), the line:
std::vector> all_boxes(tria
<https://www.dealii.org/current/doxygen/deal.II/tria__description_8cc.html#ab695339a227bd369fcea767631247cff>
.n_locally_owned_active_cells());
gets used to set the vector size. From within my subroutine, I don't have
access to the triangulation object. Am I missing a way to get
"n_locally_owned_active_cells" from the DoFHandler? Or should I just pass
this value in as an argument to my subroutine?
2) Does my approach of creating vectors to hold the bounding boxes and the
cell pointers make sense? Is there a better way to get access to the cell
once I've determined the BoundingBox index from the R-tree?
3) Somewhat related to question 1, is there a way to check that the
'closest_cell' was actually found? I had an issue with MPI where the
target point 'bc.coords' was outside the bounding box of some processes,
resulting in some errors. I set this (closest_cell_idx[0] < i-1) to check
that an actual BoundingBox was found.
Thanks for your suggestion on the R-tree. I'd welcome any critiques of my
implementation.
Cheers,
Alex
On Wednesday, November 15, 2023 at 8:35:49 AM UTC-5 Alex Quinlan wrote:
> Thanks, Wolfgang. I will abide by those guidelines.
>
> On Tuesday, November 14, 2023 at 10:15:29 PM UTC-5 Wolfgang Bangerth wrote:
>
>> On 11/14/23 16:25, Alex Quinlan wrote:
>> >
>> > I'm curious what your thoughts are on this approach. I imagine it
>> could have
>> > an advantage in some situations, depending on the size of the mesh and
>> the
>> > number of constraints to be added.
>> >
>> > I have not done any speed testing on this yet, tho. Do you think it
>> would be
>> > looking into? Or do you see some fatal flaw with this approach?
>>
>> Alex:
>> the general approach most professional programmers will ascribe to is to
>> write
>> a version of the code that is intelligible and easy to maintain. Only
>> then do
>> you worry about speed. If the code in question is fast enough (say at
>> most of
>> few percent up to 20% of the program's run time -- as measured with a
>> class
>> such as TimerOutput), then it's not worth worrying about it.
>>
>> So the questions you're asking are premature. Make the code work what it
>> is
>> you want it to do, and then you can think about its perfor
Re: [deal.II] Mixed codimensionality
Thanks, Wolfgang. Sorry for the big code dump. I've been trying to distill the problem down to find the issue. I've now eliminated the shell elements and the solid-shell coupling. In doing so, I've found a strange change to my sparsity pattern. The overlap of Process0 and Process1 is not symmetric. The PETSc error occurs in this area where I'd expect to see symmetry (at entry 36,48). [image: asym_dsp2.png] The geometry of the problem is shown in the picture below. I've colored the cells and vertices based on their ownership. [image: solid_geom.png] The dynamic sparsity pattern comes from this approach: dof_handler_sld.distribute_dofs(fe_sld); sld_owned_dofs = dof_handler_sld.locally_owned_dofs(); DynamicSparsityPattern debug_dsp(sld_owned_dofs); DoFTools::make_sparsity_pattern(dof_handler_sld, debug_dsp); It makes sense that I should have an entry at position (36,48), since dof 36 and dof 48 share a cell. Therefore, I think my stiffness matrix process is okay. My problem seems to be with my sparsity pattern. So, I'm confused as to why Proc0 would not have any Sparsity Pattern entries for (36,48) when Proc1 does indeed have (48,36). Any thoughts on where I might be causing this problem? Best regards, Alex On Thursday, February 1, 2024 at 7:56:05 PM UTC-5 Wolfgang Bangerth wrote: > On 2/1/24 14:50, Alex Quinlan wrote: > > > > Can anyone see where I may have gone wrong or what I have omitted? > > Alex: > I cannot see what the issue is, and I doubt anyone can without spending a > substantial amount of time debugging. I am slightly confused, though, by > the > error message. I *believe* that what it is saying is that you are trying > to > write into a matrix entry that is not part of the matrix -- you can only > write > into matrix entries you have previously said exist when you passed the > sparsity pattern to the matrix. So what you need to find out is who owns > row > 36, and whether the sparsity pattern you provided on that process included > entry (36,48). My best guess is that that entry was not present in the > sparsity pattern on the process that owns row 36, and then your goal needs > to > be to figure out why not -- in other words, why you did not foresee that > the > entry that you ultimately end up writing into was not included when you > constructed the sparsity pattern. > > Best > W. > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/a1cf2d0c-f86f-402d-bf3e-aeaf58758825n%40googlegroups.com.
Re: [deal.II] Mixed codimensionality
solution.reinit(dofs_per_block);
system_matrix.collect_sizes();
}
*void ElasticProblem::assemble_system()*
{
system_rhs= 0;
system_matrix = 0;
... Create FE values objects ...
for (const auto : dof_handler_sld.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_matrix_sld = 0;
cell_rhs_sld= 0;
for (const unsigned int i : fe_values_sld.dof_indices())
for (const unsigned int j : fe_values_sld.dof_indices())
for (const unsigned int q_point :
fe_values_sld.quadrature_point_indices())
{
cell_matrix_sld(i,j) = ...long equation...
}
cell->get_dof_indices(local_dof_indices_sld);
constraints.distribute_local_to_global( cell_matrix_sld, cell_rhs_sld,
local_dof_indices_sld, system_matrix.block(0,0),
system_rhs.block(0));
} // end of cell loop
MPI_Barrier(mpi_communicator);
system_matrix.compress(VectorOperation::add); *// < This is where my
issue is*
system_rhs.compress(VectorOperation::add);
// Onto Shell element assembly ...
...
...
}
The same error occurs if I skip the system_matrix.compress() at the end of
the solid elements, and instead wait until the end of the shell element
assembly.
I've looked at Step-55 and Step-70, which use the MPI::BlockSparseMatrix,
but I don't see what I'm missing. These examples use an index set
locally_relevant_dofs_per_processor, but I am using a vector of index sets,
since I have two DoF handlers.
Can anyone see where I may have gone wrong or what I have omitted?
Many thanks,
Alex
On Wednesday, January 24, 2024 at 12:00:23 PM UTC-5 Wolfgang Bangerth wrote:
>
> On 1/24/24 09:07, Alex Quinlan wrote:
> >
> > I seem to have it working now, though I've hit some other unrelated
> > snags that I need to resolve. Once I get those fixed and confirm that
> > my program is running correctly, I would be willing to work on a patch.
>
> Excellent, thank you in advance!
> Best
> W.
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/c5b355ca-72d5-4f79-9e24-011d4d28e0e1n%40googlegroups.com.
Re: [deal.II] Imported mesh problem with p4est and parallel::distributed
Thanks, Wolfgang. I'll look into it some more and eventually try with a larger problem. On Tuesday, January 30, 2024 at 11:18:01 PM UTC-5 Wolfgang Bangerth wrote: > > Alex: > > > I am running a solid-mechanics job where I import a large mesh and run > using > > parallel fullydistributed with MPI. I had been trying to run my job > using a CG > > solver with the BoomerAMG > > preconditioner (based on the example in step-40). > > > > I ran my mesh with 116,000 nodes and the solver took about 700 seconds > and > > used around 18GB of RAM. I then tried to run another version of the > same > > problem with a finer 1.5M node mesh. This one never finished because > > it's memory requirements exceeded the resources of my small cluster. > > > > Eventually, I decided to test out some different solvers and > preconditioners. > > I switched the preconditioner to Jacobi and suddenly the 116k job ran in > 40 > > seconds and only needed ~4GB of memory. > > This is too small a problem to be indicative of much. Jacobi is good for > small > problems, see also here: > > https://dealii.org/developer/doxygen/deal.II/step_6.html#step_6-Solversandpreconditioners > But it is not effective for large problems because the effort to solve > linear > problems is not O(N) with this preconditioner. AMG typically is, but it is > only better on large problems. > > I do wonder about your memory usage. A good order of magnitude for 2d > problems > is 1 kB per DoF. Perhaps twice that in 3d. You are using 30 kB. > > > > - Maybe it's because my mesh is generated externally and not refined > using the > > deal.ii mesh refinement approach? > > - Maybe it's related to the use of p:f:t instead of p:d:t? > > - Maybe I have some issue with my PETSc install? > > - Maybe I didn't properly set the parameters for the AMG preconditioner? > > All possibilities. But hard to tell on such a small problem. > > Best > W. > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/a5263d6d-b2ef-4c2b-8313-7c4bbdc5870an%40googlegroups.com.
Re: [deal.II] Imported mesh problem with p4est and parallel::distributed
Dear deal.ii community, I've had relative success running the parallel::fullydistributed, but now I've encountered some strange preconditioner/solver behavior. I am running a solid-mechanics job where I import a large mesh and run using parallel fullydistributed with MPI. I had been trying to run my job using a CG solver with the BoomerAMG preconditioner (based on the example in step-40). I ran my mesh with 116,000 nodes and the solver took about 700 seconds and used around 18GB of RAM. I then tried to run another version of the same problem with a finer 1.5M node mesh. This one never finished because it's memory requirements exceeded the resources of my small cluster. Eventually, I decided to test out some different solvers and preconditioners. I switched the preconditioner to Jacobi and suddenly the 116k job ran in 40 seconds and only needed ~4GB of memory. I was quite surprised by this. I thought that Jacobi was the bare minimum of a preconditioner and that the AMG would be better suited for this job. So, I'm curious if you have any thoughts on why my job might be reacting this way: a 20x speed-up and a 5x reduction in memory requirements. - Maybe it's because my mesh is generated externally and not refined using the deal.ii mesh refinement approach? - Maybe it's related to the use of p:f:t instead of p:d:t? - Maybe I have some issue with my PETSc install? - Maybe I didn't properly set the parameters for the AMG preconditioner? I have limited understanding of solvers and preconditioners, so I am trying to learn more about them. I'd appreciate any input you may have. Best regards, Alex On Friday, December 22, 2023 at 10:13:55 AM UTC-5 Alex Quinlan wrote: > Thanks for your response, Wolfgang. > > Your last sentence seems to be my solution. I was planning to use > parallel::fullydistributed due to the large sizes of my imported meshes > (per a suggestion on a previous post of mine: > https://groups.google.com/g/dealii/c/V5HH2pZ0Kow ) > > I wanted to run parallel:distributed primarily to compare with > parallel::fullydistributed. I was also planning to use the approach in > dealii/tests/fullydistributed_grids/copy_distributed_tria_* to convert a > p:d:t into p:d:f. > > But now, I will proceed by skipping the comparison and by using the > approach in dealii/tests/fullydistributed_grids/copy_serial_tria_*, > starting with parallel:shared:triangulation and then building the p:f:t off > of that. This allows me to side-step bug-7428 and to import a simplex mesh. > > Best regards, > Alex > > On Tuesday, December 19, 2023 at 8:26:04 PM UTC-5 Wolfgang Bangerth wrote: > >> >> Alex: >> >> You've hit on one of those bugs that every once in a while someone trips >> over, >> but that nobody really ever takes/has the time to fully debug. For sure, >> we >> would love to get some help with this issue, and the github bug report >> already >> contains a relatively small test case that should make debugging not too >> painful. The right approach is likely to use this test case, see what >> inputs >> deal.II is calling p4est's mesh creation routine, and then generate a >> stand-alone test case that only uses p4est to trigger the problem. I >> would >> expect this can be done in 200 or so lines. >> >> I believe that at the core of the problem lies a mismatch of what we >> think >> p4est wants, and what p4est really expects. This also illustrates the >> problem >> with debugging the issue: Those of us who know deal.II don't know p4est, >> and >> the other way around. It requires someone willing to create such a >> stand-alone >> test case, and then talk to the p4est folks about why this is wrong. Once >> we >> understand why the input is wrong/unacceptable to p4est, it should be >> relatively straightforward to come up with a way to give p4est what it >> wants. >> The *understanding* of the problem is the key step. >> >> >> > Can you help me out by answering a couple questions? >> > >> > * Is there a work-around for this issue? >> >> No. >> >> > o Is there an alternative to p4est? >> >> In the long run, the p4est people are working on a library called tetcode >> that >> extends p4est. But I don't think any of us have ever tried to see what it >> would take to replace p4est by tetcode. >> >> >> > o Are there modifications I can make to my mesh/program to avoid this >> issue? >> >> Not sure. Absent an understanding of why p4est is upset, any >> modifications >> would just be poking in the dark. >> >> >> > * What would be involved
Re: [deal.II] Mixed codimensionality
Thanks, Wolfgang. I seem to have it working now, though I've hit some other unrelated snags that I need to resolve. Once I get those fixed and confirm that my program is running correctly, I would be willing to work on a patch. Thanks, Alex On Friday, January 19, 2024 at 1:00:51 PM UTC-5 Wolfgang Bangerth wrote: > > On 1/18/24 12:56, Alex Quinlan wrote: > > Where I start to get confused is when I try to reinitialize the system > > matrix. What I want to do is something like: > > > > system_matrix.reinit(dsp, mpi_communicator); > > > > > > This has my MPI communicator and the sparsity pattern that I’ve built > > up. However, this isn’t a valid call to > > PETScWrappers::MPI::BlockSparseMatrix.reinit(). There’s a similar > > function that takes the arguments ( *const std::vector< IndexSet > > <https://www.dealii.org/current/doxygen/deal.II/classIndexSet.html> > > ,*const BlockDynamicSparsityPattern < > https://www.dealii.org/current/doxygen/deal.II/classBlockDynamicSparsityPattern.html> > > , const MPI_Comm < > https://www.dealii.org/current/doxygen/deal.II/classMPI__Comm.html> com ) > > > > > > I don’t really understand what I would put in for the "sizes" vector. > > What exactly am I trying to pass with this argument? Is it all of the > > locally owned/relevant dofs? Do I just combine the vector of locally > > owned shell dofs and locally owned solid dofs? > > The argument is poorly named. It is a vector of index sets (in your case > a vector of 2 index sets) each of which contains which of the elements > of a block are locally owned. Perhaps the corresponding function in > other block classes uses a better named and/or documented argument. (The > name dates back to a time when we just passed down how many rows each > process owns, rather than being explicit -- with an index set -- which > rows these actually are.) > > This would probably be worth fixing. Would you be willing to write a > patch that renames the argument? > > Best > W. > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/3d597190-6c08-4250-a3e6-a895aa4aa4dbn%40googlegroups.com.
Re: [deal.II] Mixed codimensionality
Thanks for you response, Wolfgang. The triangulation partitioning makes a
lot of sense now.
In that case, perhaps my problem is not too difficult. I am getting caught
up with the PETScWrappers::MPI::BlockSparseMatrix.reinit() method.
This is my approach:
Distributing dofs and collecting the locally owned and locally relevant
sets.
// Solid
dof_handler_sld.distribute_dofs(fe_sld);
sld_owned_dofs = dof_handler_sld.locally_owned_dofs();
sld_relevant_dofs =
DoFTools::extract_locally_relevant_dofs(dof_handler_sld);
// Shell
dof_handler_sh.distribute_dofs(fe_sh);
sh_owned_dofs = dof_handler_sh.locally_owned_dofs();
sh_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler_sh);
Then I apply boundary conditions to both solid and shell dof_handlers. Then
I define the coupling points and retrieve the coupled DOFs:
sld_coup_dofs = nodal_coupling(coup_points , dof_handler_sld,
tria_pft_sld.n_vertices() );
sh_coup_dofs = nodal_coupling(coup_points , dof_handler_sh,
tria_pft_sh.n_vertices() );
const std::vector dofs_per_block = {n_dof_sld, n_dof_sh};
const std::vector locally_owned_sizes =
{sld_owned_dofs.size(), sh_owned_dofs.size()};
I reinitialize each of the blocks in my sparsity pattern:
BlockDynamicSparsityPattern dsp(dofs_per_block , dofs_per_block);
for (unsigned int i=0; ihttps://www.dealii.org/current/doxygen/deal.II/classIndexSet.html> >
,* const BlockDynamicSparsityPattern
<https://www.dealii.org/current/doxygen/deal.II/classBlockDynamicSparsityPattern.html>
, const MPI_Comm
<https://www.dealii.org/current/doxygen/deal.II/classMPI__Comm.html> com )
I don’t really understand what I would put in for the "sizes" vector. What
exactly am I trying to pass with this argument? Is it all of the locally
owned/relevant dofs? Do I just combine the vector of locally owned shell
dofs and locally owned solid dofs?
Thanks,
Alex
On Wednesday, January 17, 2024 at 8:24:44 PM UTC-5 Wolfgang Bangerth wrote:
> On 1/17/24 14:40, Alex Quinlan wrote:
> >
> > If I do this, then it seems that I'd end up with two conflicting
> > partitions. If I have n_mpi_processes=8 , for example, this approach
> > would assign ~1/8th of the solid mesh to process 1 AND ~1/8th of the
> > shell mesh to process 1. This seems like a problem.
>
> No, this is exactly what you want. That's because you will have phases
> in your program when you work on the solid mesh, and you want all 8
> processes to participate in that. And then you have a phase where you
> going to do something on the surface mesh and, again, you want to have
> all 8 processes work on that.
>
> It is *so* much harder to write software in which different processes do
> different things at the same time. Don't go there.
>
>
> > As you pointed out, I will have an issue with the off-diagonal matrices
> > and the coupling between the two triangulations. Is it possible for me
> > to build up all of the sparsity patterns and matrices /and then /perform
> > the partitioning?
>
> You could of course, but then what's the point of using multiple
> processors if you do the bulk of the work on only one process?
>
> Best
> W.
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/6c55938a-2132-4747-ac3d-712669ded3a1n%40googlegroups.com.
Re: [deal.II] Mixed codimensionality
Hi Wolfgang,
Thanks for your feedback. I have this two-triangulation problem working in
serial, and I have other jobs that work with MPI. When I put them
together, I end up with something like the code below. It seems like there
will be some issues, though.
--- Reading in two triangulations
std::ifstream in;
in.open("input/solid.inp");
GridIn<3> grid_in_sld;
grid_in_sld.attach_triangulation(triangulation_sld);
grid_in_sld.read_abaqus(in);
std::ifstream in2;
in2.open("input/flat3x3.inp");
GridIn<2,3> grid_in_sh;
grid_in_sh.attach_triangulation(triangulation_sh);
grid_in_sh.read_abaqus(in2);
--- Partitioning the two meshes
GridTools::partition_triangulation_zorder(
Utilities::MPI::n_mpi_processes(mpi_communicator),
triangulation_sld); // maybe divide n_mpi by 2
GridTools::partition_multigrid_levels(triangulation_sld);
auto construction_data_sld =
TriangulationDescription::Utilities::create_description_from_triangulation(
triangulation_sld, mpi_communicator,
TriangulationDescription::Settings::construct_multigrid_hierarchy);
tria_pft_sld.create_triangulation(construction_data_sld);
GridTools::partition_triangulation_zorder(
Utilities::MPI::n_mpi_processes(mpi_communicator),
triangulation_sh);// maybe divide n_mpi by 2
GridTools::partition_multigrid_levels(triangulation_sh);
auto construction_data_sh =
TriangulationDescription::Utilities::create_description_from_triangulation(
triangulation_sh, mpi_communicator,
TriangulationDescription::Settings::construct_multigrid_hierarchy);
tria_pft_sh.create_triangulation(construction_data_sh);
If I do this, then it seems that I'd end up with two conflicting
partitions. If I have n_mpi_processes=8 , for example, this approach would
assign ~1/8th of the solid mesh to process 1 AND ~1/8th of the shell mesh
to process 1. This seems like a problem.
I was thinking that maybe I could set it up such that the larger solid mesh
could be split into 6 partitions over processes 0-5. Then the smaller
shell mesh could be split into 2 partitions using processes 6 and 7.
However, I didn't see any way to have the process numbers start at a
non-zero number.
As you pointed out, I will have an issue with the off-diagonal matrices and
the coupling between the two triangulations. Is it possible for me to
build up all of the sparsity patterns and matrices *and then *perform the
partitioning?
Thanks,
Alex
On Tuesday, January 16, 2024 at 10:04:52 PM UTC-5 Wolfgang Bangerth wrote:
> On 1/16/24 07:03, Alex Quinlan wrote:
> >
> > I've come up against another issue in my mixed-dimensional mesh. I am
> now
> > trying to implement MPI, which requires partitioning the triangulation.
> >
> > In my case, I have two triangulations: one for the solid <3,3> elements
> and
> > one for the shell <2,3> elements. These two triangulations are couple
> at
> > multiple points. I combine them into a blockmatrix with the form:
> >
> > [ A_solid ,B_coupling;
> > C_coupling,D_shell ]
> >
> > I may also need to expand this in the future to a matrix with 3x3 or nxn
> blocks.
> >
> > Do you have any thoughts on how I can partition this problem? The suite
> of
> > GridTools::partition_triangulation seem like they will not work for my
> > two-triangulation approach.
>
> Can you elaborate what the problem you are facing is? At least for the
> diagonal blocks of the matrix, you would simply use a matrix class such as
> PETScWrappers::MPI::SparseMatrix that knows how to partition itself onto
> different processes. Similarly, for your vectors, you'd use the
> appropriate
> block vector class in which each vector block is partitioned. You'd have
> to
> initialize these individually, but I suspect you're already doing that
> anyway.
> In each case, blocking (into the 2x2 structure you show above) and
> partitioning (among MPI processes) are orthogonal issues.
>
> The biggest issue I see is how you can build the off-diagonal blocks above
> because to do so, you need to know information about both meshes on the
> same
> process, and generally the partition of one triangulation owned by a
> specific
> process does not geometrically overlap with the partition of the other
> triangulation owned by the same process.
>
> 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
htt
Re: [deal.II] Mixed codimensionality
Dear all, I've come up against another issue in my mixed-dimensional mesh. I am now trying to implement MPI, which requires partitioning the triangulation. In my case, I have two triangulations: one for the solid <3,3> elements and one for the shell <2,3> elements. These two triangulations are couple at multiple points. I combine them into a blockmatrix with the form: [ A_solid ,B_coupling; C_coupling,D_shell ] I may also need to expand this in the future to a matrix with 3x3 or nxn blocks. Do you have any thoughts on how I can partition this problem? The suite of GridTools::partition_triangulation seem like they will not work for my two-triangulation approach. Best regards, Alex On Thursday, November 9, 2023 at 1:24:29 PM UTC-5 Alex Quinlan wrote: > Sorry, please disregard my last question. I've reviewed the Step-20 solver > that I tried to use and I see that it is a different set-up to my problem. > I'm now working to apply the Schur complement in the correct manner. > > On Wednesday, November 8, 2023 at 10:32:35 AM UTC-5 Alex Quinlan wrote: > >> Alright, I've had some decent success, though I've hit another snag in >> the solver section. Your suggestions have opened the doors and actually >> allowed me to shift gears to my preferred approach. >> >> I've changed the end of Setup_System() to be: >> >> dsp.collect_sizes(); >> >> *DoFTools::make_sparsity_pattern(dof_handler_sld, dsp.block(0,0)); >> DoFTools::make_sparsity_pattern(dof_handler_sh, dsp.block(1,1));* >> >> sparsity_pattern.copy_from(dsp); >> system_matrix.reinit(sparsity_pattern); >> system_rhs.reinit(dofs_per_block); >> solution.reinit(dofs_per_block); >> std::cout << " " << std::endl; >> >> Doing so allowed me to successfully run your suggested code: >> >> for (unsigned int i = 0; i < dofs_per_cell_sld; ++i) >> for (unsigned int j = 0; j < dofs_per_cell_sld; ++j) >> system_matrix.*block(0,0)*.add(local_dof_indices_sld[i], >> >> local_dof_indices_sld[j], >> cell_matrix_sld(i, j)); >> >> >> This success prompted me to try a different approach to applying the cell >> stiffness matrix to the global stiffness matrix. I have an object of >> AffineConstraints, so I wanted to use the >> AffineConstraints::distribute_local_to_global(). >> >> For the block system, I prepared the constraints like this: (Note: >> nodal_bcs() is a routine that I've written to apply constraints that I've >> read in from a .csv file) >> >> // Set up Solid Constraints >> constraints_sld.clear(); >> nodal_bcs<3,3>(solidbcfile, dof_handler_sld ,constraints_sld); >> constraints_sld.close(); >> >> // Set-up Shell Constraints >> constraints_sh.clear(); >> nodal_bcs<2,3>(shellbcfile, dof_handler_sh ,constraints_sh); >> constraints_sh.close(); >> >> // Combine Solid and Shell constraints >> constraints.copy_from(constraints_sld); >> constraints_sh.shift(n_dof_sld); >> constraints.merge(constraints_sh ); >> >> This allowed me to add this line to the Cell loop in the >> Assemble_System(): >> >> constraints.distribute_local_to_global( >> cell_matrix_sld, cell_rhs_sld, local_dof_indices_sld, >> system_matrix.block(0,0), system_rhs.block(0)); >> >> >> All of this now seems to work and the matrices seem to be correct. >> >> - >> >> Now, I'm having trouble with the solver. I think the cause is the empty >> off-diagonal matrix I am using the approach outlined in Step-20: >> >>const auto = system_matrix.block(0, 0); >> const auto = system_matrix.block(0, 1); >> const auto = system_rhs.block(0); >> const auto = system_rhs.block(1); >> auto = solution.block(0); >> auto = solution.block(1); >> const auto op_M = linear_operator(M); >> const auto op_B = linear_operator(B); >> >> ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10); >> SolverCG> solver_M(reduction_control_M); >> PreconditionJacobi> preconditioner_M; >> preconditioner_M.initialize(M); >> >> const auto op_M_inv = inverse_operator(op_M, solver_M, >> preconditioner_M); >> const auto op_S = transpose_operator(op_B) * op_M_inv * op_B; >> const auto op_aS = transpose_operator(op_B) * >> linear_operato
Re: [deal.II] Imported mesh problem with p4est and parallel::distributed
Thanks for your response, Wolfgang. Your last sentence seems to be my solution. I was planning to use parallel::fullydistributed due to the large sizes of my imported meshes (per a suggestion on a previous post of mine: https://groups.google.com/g/dealii/c/V5HH2pZ0Kow ) I wanted to run parallel:distributed primarily to compare with parallel::fullydistributed. I was also planning to use the approach in dealii/tests/fullydistributed_grids/copy_distributed_tria_* to convert a p:d:t into p:d:f. But now, I will proceed by skipping the comparison and by using the approach in dealii/tests/fullydistributed_grids/copy_serial_tria_*, starting with parallel:shared:triangulation and then building the p:f:t off of that. This allows me to side-step bug-7428 and to import a simplex mesh. Best regards, Alex On Tuesday, December 19, 2023 at 8:26:04 PM UTC-5 Wolfgang Bangerth wrote: > > Alex: > > You've hit on one of those bugs that every once in a while someone trips > over, > but that nobody really ever takes/has the time to fully debug. For sure, > we > would love to get some help with this issue, and the github bug report > already > contains a relatively small test case that should make debugging not too > painful. The right approach is likely to use this test case, see what > inputs > deal.II is calling p4est's mesh creation routine, and then generate a > stand-alone test case that only uses p4est to trigger the problem. I would > expect this can be done in 200 or so lines. > > I believe that at the core of the problem lies a mismatch of what we think > p4est wants, and what p4est really expects. This also illustrates the > problem > with debugging the issue: Those of us who know deal.II don't know p4est, > and > the other way around. It requires someone willing to create such a > stand-alone > test case, and then talk to the p4est folks about why this is wrong. Once > we > understand why the input is wrong/unacceptable to p4est, it should be > relatively straightforward to come up with a way to give p4est what it > wants. > The *understanding* of the problem is the key step. > > > > Can you help me out by answering a couple questions? > > > > * Is there a work-around for this issue? > > No. > > > o Is there an alternative to p4est? > > In the long run, the p4est people are working on a library called tetcode > that > extends p4est. But I don't think any of us have ever tried to see what it > would take to replace p4est by tetcode. > > > > o Are there modifications I can make to my mesh/program to avoid this > issue? > > Not sure. Absent an understanding of why p4est is upset, any modifications > would just be poking in the dark. > > > > * What would be involved with fixing the bug? > > See above. > > > > On a possibly related issue: > > > > * What is needed to allow for simplex meshes with > parallel::distributed? Is > > this also a connectivity issue? > > p::d::Triangulation builds on p4est, which is exclusively for quad/hex > meshes. > Switching to tetcode would allow for the use of other cell types as well, > but > I don't think any of us have given substantial thought to what such a > switch > would require. > > There is the parallel::fullydistributed::Triangulation class, which I > *think* > works with simplex meshes. It is not covered by tutorial programs yet, > though. > > Best > W. > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/4ef98ac2-983c-4181-b930-54f6e4912ad5n%40googlegroups.com.
[deal.II] Imported mesh problem with p4est and parallel::distributed
Dear deal.ii community, I am working with a mesh that is imported via a modified version of GridIn::read_abaqus(). I'm able to run my mesh and job with parallel::shared without any issues. However, when I go to use parallel::distributed, I run into an issue with p8est connectivity: void dealii::parallel::distributed::Triangulation::copy_new_triangulation_to_p4est(std::integral_constant) [with int dim = 3; int spacedim = 3] The violated condition was: p8est_connectivity_is_valid(connectivity) == 1 It seems that this is a known issue: https://github.com/dealii/dealii/issues/7428 I've read through the issue thread, but I don't fully understand what is going on. It's not clear to me if this is a bug with p4est or with deal.ii. There is a listed work-around that seems to rearrange the nodes, but I'm not sure if that's applicable in my situation. Can you help me out by answering a couple questions? - Is there a work-around for this issue? - Is there an alternative to p4est? - Are there modifications I can make to my mesh/program to avoid this issue? - What would be involved with fixing the bug? On a possibly related issue: - What is needed to allow for simplex meshes with parallel::distributed? Is this also a connectivity issue? Many thanks, Alex -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/d8ee7ec8-4ca9-4c3a-82f6-33b7c356cd9an%40googlegroups.com.
Re: [deal.II] Artificial cells on large imported meshes during parallel::distributed
Thanks, Daniel. I will look into parallel::fullydistributed. It appears that there are no step-## tutorials with parallel::fully distributed, and there's less documentation compared to parallel::distributed. I found this repo that looks like it has some good examples about parallel::fullydistributed, so I will start here: https://github.com/peterrum/dealii-pft Best regards, Alex On Wednesday, December 13, 2023 at 11:27:41 AM UTC-5 d.arnd...@gmail.com wrote: > Alex, > > in this case, the triangulation is copied for all processes when using > parallel::distributed::Triangulation. You might want to look into using > parallel::fullydistributed::Triangulation if your unrefined mesh is very > fine. > > Best, > Daniel > > On Wed, Dec 13, 2023 at 8:58 AM Alex Quinlan wrote: > >> Dear Deal.ii Community, >> >> I have been reading the "Parallel computing with multiple processors >> using distributed memory >> <https://www.dealii.org/current/doxygen/deal.II/group__distributed.html>" >> module and I'd like to clarify something about the artificial cells shown >> in the picture below. >> >> [image: parallel_dist.png] >> >> I understand from the module that each processor will hold onto it's >> local portion of cells, along with the "ghost cells" on the border with >> other processors. It also includes the dark blue *artificial cells*, >> which "ensure that each processor has a mesh that has all the coarse level >> cells and that respects the invariant that neighboring cells can not differ >> by more than one level of refinement." >> >> I believe I understand how this works for a problem that begins with a >> coarse mesh and then undergoes a series of mesh refinements. In this case, >> a processor stores it's local refined mesh plus a coarser mesh for >> non-relevant cells. >> >> However, I'm unclear about how this works if I were to import a fully >> refined mesh. Each processor would have it's local mesh, but would it also >> be storing the artificial cells at the same refinement level as the >> imported mesh? Or is there some method of coarsening the mesh on the >> artificial cells? >> >> So, for example, if I were to import this mesh:[image: meshfine.png] >> >> will my processor be stuck with a mesh of artificial cells like this? >> >> CASE 1) >> [image: fine-partition.png] >> >> Or is there a way that the mesh will auto-coarsen the artificial cells to >> something preferred like this? >> >> CASE 2) >> [image: coarse-partition.png] >> (note: I realize that the ghost cells are not represented properly in >> case 2, and that they should be the same as case 1. ) >> >> >> I am concerned that CASE 1 will cause memory problems when I have very >> large imported meshes. Can you enlighten me on this issue? >> >> Best regards, >> Alex >> >> -- >> The deal.II project is located at http://www.dealii.org/ >> For mailing list/forum options, see >> https://groups.google.com/d/forum/dealii?hl=en >> --- >> You received this message because you are subscribed to the Google Groups >> "deal.II User Group" group. >> To unsubscribe from this group and stop receiving emails from it, send an >> email to dealii+un...@googlegroups.com. >> To view this discussion on the web visit >> https://groups.google.com/d/msgid/dealii/b7dc13e7-7b99-4be3-8e23-d1450862cecfn%40googlegroups.com >> >> <https://groups.google.com/d/msgid/dealii/b7dc13e7-7b99-4be3-8e23-d1450862cecfn%40googlegroups.com?utm_medium=email_source=footer> >> . >> > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/8648882f-fc10-4968-8f98-207ad3d324bcn%40googlegroups.com.
[deal.II] Reporting issue with repo
Dear deal.ii devs, I cloned the deal.ii repo yesterday and found that the header file grid_reorder.h has been removed from include/deal.II/grid/. I searched the logs and the documentation and it appeared that the file *should* still be there, though I'm not 100% sure of that. I am not very savvy with git, so I thought I would just post here. Is this the correct place to make this kind of report, or is there a preferred process that I should take in the future? Thanks, Alex -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/3c1a7409-a834-4fed-b1b5-48caa5de7a69n%40googlegroups.com.
Re: [deal.II] Importing nodal BCs and accessing vertices
Thanks, Wolfgang. I will abide by those guidelines. On Tuesday, November 14, 2023 at 10:15:29 PM UTC-5 Wolfgang Bangerth wrote: > On 11/14/23 16:25, Alex Quinlan wrote: > > > > I'm curious what your thoughts are on this approach. I imagine it could > have > > an advantage in some situations, depending on the size of the mesh and > the > > number of constraints to be added. > > > > I have not done any speed testing on this yet, tho. Do you think it > would be > > looking into? Or do you see some fatal flaw with this approach? > > Alex: > the general approach most professional programmers will ascribe to is to > write > a version of the code that is intelligible and easy to maintain. Only then > do > you worry about speed. If the code in question is fast enough (say at most > of > few percent up to 20% of the program's run time -- as measured with a > class > such as TimerOutput), then it's not worth worrying about it. > > So the questions you're asking are premature. Make the code work what it > is > you want it to do, and then you can think about its performance. > > Best > W. > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/08be57a0-9ea9-4726-96e1-22087ca5c3e7n%40googlegroups.com.
Re: [deal.II] Importing nodal BCs and accessing vertices
Hi Wolfgang, I've come back to this question after discovering some Grid tools. I am now assessing GridTools::find_closest_vertex and GridTools::find_active_cell_around_point const Point<3> pt1(0.0, 15.0, 8.20754); std::vector marked_vertices(triangulation.n_vertices(), true); // Find the vertex closest to the point unsigned int closest_v1 = GridTools::find_closest_vertex( dof_handler, pt1); // Mark the vertex for a faster cell search marked_vertices[closest_v1] = true; // Find the cell containing the point (or a cell containing, if there are multiple) auto cell_and_ref_point = GridTools::find_active_cell_around_point( dof_handler,pt1, marked_vertices, 1.e-4); // Search the vertices, then apply constraints for (unsigned int v = 0; vn_vertices() ; v++) // might be able to use 'vertex_iterator' if (cell_and_ref_point->vertex(v).distance(load1) < 1.e-4) constraints.add_line(cell_and_ref_point->vertex_dof_index(v,0, cell_and_ref_point->active_fe_index() )); I'm curious what your thoughts are on this approach. I imagine it could have an advantage in some situations, depending on the size of the mesh and the number of constraints to be added. I have not done any speed testing on this yet, tho. Do you think it would be looking into? Or do you see some fatal flaw with this approach? Thanks, Alex On Wednesday, January 18, 2023 at 8:31:11 AM UTC-5 Alex Quinlan wrote: > Thanks Wolfgang. I appreciate the feedback. I'll see if I can implement > some of the cost-saving tips that you suggested. > > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/ce083242-20b0-4091-b93a-142ddcc2534bn%40googlegroups.com.
Re: [deal.II] Mixed codimensionality
Sorry, please disregard my last question. I've reviewed the Step-20 solver that I tried to use and I see that it is a different set-up to my problem. I'm now working to apply the Schur complement in the correct manner. On Wednesday, November 8, 2023 at 10:32:35 AM UTC-5 Alex Quinlan wrote: > Alright, I've had some decent success, though I've hit another snag in the > solver section. Your suggestions have opened the doors and actually > allowed me to shift gears to my preferred approach. > > I've changed the end of Setup_System() to be: > > dsp.collect_sizes(); > > *DoFTools::make_sparsity_pattern(dof_handler_sld, dsp.block(0,0)); > DoFTools::make_sparsity_pattern(dof_handler_sh, dsp.block(1,1));* > > sparsity_pattern.copy_from(dsp); > system_matrix.reinit(sparsity_pattern); > system_rhs.reinit(dofs_per_block); > solution.reinit(dofs_per_block); > std::cout << " " << std::endl; > > Doing so allowed me to successfully run your suggested code: > > for (unsigned int i = 0; i < dofs_per_cell_sld; ++i) > for (unsigned int j = 0; j < dofs_per_cell_sld; ++j) > system_matrix.*block(0,0)*.add(local_dof_indices_sld[i], > > local_dof_indices_sld[j], > cell_matrix_sld(i, j)); > > > This success prompted me to try a different approach to applying the cell > stiffness matrix to the global stiffness matrix. I have an object of > AffineConstraints, so I wanted to use the > AffineConstraints::distribute_local_to_global(). > > For the block system, I prepared the constraints like this: (Note: > nodal_bcs() is a routine that I've written to apply constraints that I've > read in from a .csv file) > > // Set up Solid Constraints > constraints_sld.clear(); > nodal_bcs<3,3>(solidbcfile, dof_handler_sld ,constraints_sld); > constraints_sld.close(); > > // Set-up Shell Constraints > constraints_sh.clear(); > nodal_bcs<2,3>(shellbcfile, dof_handler_sh ,constraints_sh); > constraints_sh.close(); > > // Combine Solid and Shell constraints > constraints.copy_from(constraints_sld); > constraints_sh.shift(n_dof_sld); > constraints.merge(constraints_sh ); > > This allowed me to add this line to the Cell loop in the Assemble_System(): > > constraints.distribute_local_to_global( > cell_matrix_sld, cell_rhs_sld, local_dof_indices_sld, > system_matrix.block(0,0), system_rhs.block(0)); > > > All of this now seems to work and the matrices seem to be correct. > > - > > Now, I'm having trouble with the solver. I think the cause is the empty > off-diagonal matrix I am using the approach outlined in Step-20: > >const auto = system_matrix.block(0, 0); > const auto = system_matrix.block(0, 1); > const auto = system_rhs.block(0); > const auto = system_rhs.block(1); > auto = solution.block(0); > auto = solution.block(1); > const auto op_M = linear_operator(M); > const auto op_B = linear_operator(B); > > ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10); > SolverCG> solver_M(reduction_control_M); > PreconditionJacobi> preconditioner_M; > preconditioner_M.initialize(M); > > const auto op_M_inv = inverse_operator(op_M, solver_M, > preconditioner_M); > const auto op_S = transpose_operator(op_B) * op_M_inv * op_B; > const auto op_aS = transpose_operator(op_B) * > linear_operator(preconditioner_M) * op_B; > > IterationNumberControl iteration_number_control_aS(30, 1.e-18); > SolverCG> solver_aS(iteration_number_control_aS); > const auto preconditioner_S = inverse_operator(op_aS, solver_aS, > PreconditionIdentity()); > const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G; > SolverControlsolver_control_S(2000, 1.e-12); > SolverCG> solver_S(solver_control_S); > const auto op_S_inv = inverse_operator(op_S, solver_S, > preconditioner_S); > > * P = op_S_inv * schur_rhs;* > std::cout << solver_control_S.last_step() > << " CG Schur complement iterations to obtain convergence." > << std::endl; > U = op_M_inv * (F - op_B * P); > > At the red line above, I am given this error: > > An error occurred in line <637> of file in > function > void dealii::internal::SolverCG::IterationWorker MatrixType, PreconditionerType, > >::do_iteration(unsigned int) [with VectorType = dealii::Vector; > MatrixType = dealii::LinearOperator, > dealii::Vector, > de
Re: [deal.II] Mixed codimensionality
Alright, I've had some decent success, though I've hit another snag in the solver section. Your suggestions have opened the doors and actually allowed me to shift gears to my preferred approach. I've changed the end of Setup_System() to be: dsp.collect_sizes(); *DoFTools::make_sparsity_pattern(dof_handler_sld, dsp.block(0,0)); DoFTools::make_sparsity_pattern(dof_handler_sh, dsp.block(1,1));* sparsity_pattern.copy_from(dsp); system_matrix.reinit(sparsity_pattern); system_rhs.reinit(dofs_per_block); solution.reinit(dofs_per_block); std::cout << " " << std::endl; Doing so allowed me to successfully run your suggested code: for (unsigned int i = 0; i < dofs_per_cell_sld; ++i) for (unsigned int j = 0; j < dofs_per_cell_sld; ++j) system_matrix.*block(0,0)*.add(local_dof_indices_sld[i], local_dof_indices_sld[j], cell_matrix_sld(i, j)); This success prompted me to try a different approach to applying the cell stiffness matrix to the global stiffness matrix. I have an object of AffineConstraints, so I wanted to use the AffineConstraints::distribute_local_to_global(). For the block system, I prepared the constraints like this: (Note: nodal_bcs() is a routine that I've written to apply constraints that I've read in from a .csv file) // Set up Solid Constraints constraints_sld.clear(); nodal_bcs<3,3>(solidbcfile, dof_handler_sld ,constraints_sld); constraints_sld.close(); // Set-up Shell Constraints constraints_sh.clear(); nodal_bcs<2,3>(shellbcfile, dof_handler_sh ,constraints_sh); constraints_sh.close(); // Combine Solid and Shell constraints constraints.copy_from(constraints_sld); constraints_sh.shift(n_dof_sld); constraints.merge(constraints_sh ); This allowed me to add this line to the Cell loop in the Assemble_System(): constraints.distribute_local_to_global( cell_matrix_sld, cell_rhs_sld, local_dof_indices_sld, system_matrix.block(0,0), system_rhs.block(0)); All of this now seems to work and the matrices seem to be correct. - Now, I'm having trouble with the solver. I think the cause is the empty off-diagonal matrix I am using the approach outlined in Step-20: const auto = system_matrix.block(0, 0); const auto = system_matrix.block(0, 1); const auto = system_rhs.block(0); const auto = system_rhs.block(1); auto = solution.block(0); auto = solution.block(1); const auto op_M = linear_operator(M); const auto op_B = linear_operator(B); ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10); SolverCG> solver_M(reduction_control_M); PreconditionJacobi> preconditioner_M; preconditioner_M.initialize(M); const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M); const auto op_S = transpose_operator(op_B) * op_M_inv * op_B; const auto op_aS = transpose_operator(op_B) * linear_operator(preconditioner_M) * op_B; IterationNumberControl iteration_number_control_aS(30, 1.e-18); SolverCG> solver_aS(iteration_number_control_aS); const auto preconditioner_S = inverse_operator(op_aS, solver_aS, PreconditionIdentity()); const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G; SolverControlsolver_control_S(2000, 1.e-12); SolverCG> solver_S(solver_control_S); const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S); * P = op_S_inv * schur_rhs;* std::cout << solver_control_S.last_step() << " CG Schur complement iterations to obtain convergence." << std::endl; U = op_M_inv * (F - op_B * P); At the red line above, I am given this error: An error occurred in line <637> of file in function void dealii::internal::SolverCG::IterationWorker >::do_iteration(unsigned int) [with VectorType = dealii::Vector; MatrixType = dealii::LinearOperator, dealii::Vector, dealii::internal::LinearOperatorImplementation::EmptyPayload>; PreconditionerType = dealii::LinearOperator, dealii::Vector, dealii::internal::LinearOperatorImplementation::EmptyPayload>; = int] The violated condition was: std::abs(p_dot_A_dot_p) != 0. Additional information: A piece of code is attempting a division by zero. This is likely going to lead to results that make no sense. I'm guessing this is related to the inverse operators? So, I'm curious if there's some way to use this approach for problems that have a Schur complements that can be either zero or non-zero. Thanks, Alex On Wednesday, November 8, 2023 at 9:50:17 AM UTC-5 Alex Quinlan wrote: > Dear Luca and Wolfgang, > > Thanks for your pointers; I will work to implement them. > > My github account is: https://github.com/AlexQuinlan >
Re: [deal.II] Mixed codimensionality
Dear Luca and Wolfgang,
Thanks for your pointers; I will work to implement them.
My github account is: https://github.com/AlexQuinlan
Thanks,
Alex
On Wednesday, November 8, 2023 at 3:36:43 AM UTC-5 luca@gmail.com wrote:
> Dear Alex,
>
> if you send me your github user, I’ll give you access to a code that does
> exactly what you are trying to do, so that you can get some inspiration on
> how to assemble the coupling matrices.
>
> There is an open PR (the review has been stalling for some time on that)
> https://github.com/dealii/dealii/pull/15773 that tackle exactly this use
> case.
>
> Let me know if you need some assistance.
>
> L.
>
> > On 8 Nov 2023, at 05:02, Wolfgang Bangerth
> wrote:
> >
> >
> > Alex:
> >
> > > [...]
> >
> > I think everything up to this point looks reasonable.
> >
> >> From what I can see, the system_matrix looks right with n_blocks = 2
> and start_indices = {0, 24, 72}, as expected.
> >> Next, I run through the Assemble_System() to create a local cell
> matrix. I start with the Solid elements and cell_matrix_sld. I'm skipping
> the calculation of the cell matrix for brevity's sake.
> >> for (const auto : dof_handler_sld.active_cell_iterators())
> >> {
> >> // *** Placeholder for calculating cell_matrix_sld, a 24x24 matrix ***
> >> cell->get_dof_indices(local_dof_indices_sld);
> >> // Add the cell matrix to the system matrix
> >> for (unsigned int i = 0; i < dofs_per_cell_sld; ++i)
> >> for (unsigned int j = 0; j < dofs_per_cell_sld; ++j)
> >> system_matrix.add(local_dof_indices_sld[i],
> >> local_dof_indices_sld[j],
> >> cell_matrix_sld(i, j));
> >> } // end of cell loop
> >> The bit of code 'system_matrix.add(...)' works for when i=0 and j=0,
> but fails on the next loop when j=1. I get this error:
> >
> > The issue here is that you're indexing into the *global* matrix with the
> indices of only one of the DoFHandler objects. That can't go well. It
> *should* work if you do the indexing only in the matrix block that
> corresponds to dof_handler_sld:
> >
> > system_matrix.block(block_index_sld, block_index_sld)
> > .add(local_dof_indices_sld[i],
> > local_dof_indices_sld[j],
> > cell_matrix_sld(i, j));
> >
> > where you define block_index_sld somewhere to either be zero or one.
> >
> > In your case, this won't work yet because you hadn't built the sparsity
> pattern yet. You can do that the same way:
> >
> > DoFTools::make_sparsity_pattern(dof_handler_sld,
> > block_dsp.block(block_index_sld, block_index_sld));
> >
> > Best
> > W.
> >
> > --
> >
> > Wolfgang Bangerth email: bang...@colostate.edu
> > www: http://www.math.colostate.edu/~bangerth/
> >
> >
> > --
> > The deal.II project is located at http://www.dealii.org/
> > For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> > --- You received this message because you are subscribed to the Google
> Groups "deal.II User Group" group.
> > To unsubscribe from this group and stop receiving emails from it, send
> an email to dealii+un...@googlegroups.com.
> > To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/9a890d55-91f6-ec16-98be-acdb730919d3%40colostate.edu
> .
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/f30c3007-e56d-4beb-9a0c-1a1b6f77474bn%40googlegroups.com.
Re: [deal.II] Mixed codimensionality
Dear Wolfgang and Daniel,
Thanks for your directions. I've reviewed Step-20, Step-60, and Step-70
and started working on my own example.
In this case, I'm trying to combine two types of dimensionality: the <3,3>
Solid Elements and the <2,3> "Shell" Elements. I've attached a picture of
the meshes. I am trying to do a block system. For the moment, I am not
trying to link the two meshes (off-diagonal blocks), only get them into the
same system matrix. Each is independently constrained with appropriate
boundary conditions. I've run the meshes separately and they work fine.
My understanding is that I cannot combine several objects because of the
difference in . So, I've created duplicates of several
objects, using "_sld" to note for the Solid and "_sh" for the shell.
Triangulation<3> triangulation_sld;
Triangulation<2,3> triangulation_sh;
DoFHandler<3>dof_handler_sld;
DoFHandler<2,3>dof_handler_sh;
FESystem<3> fe_sld;
FESystem<2,3> fe_sh;
const MappingFE<3>mapping_sld;
const MappingFE<2,3> mapping_sh;
const QGauss<3> quadrature_sld;
const QGauss<2> quadrature_sh;
I then create the following objects that will be shared by both the Solid
and Shell elements.
AffineConstraints constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix system_matrix;
BlockSparseMatrix preconditioner_matrix;
BlockVector solution;
BlockVector system_rhs;
BlockVector locally_relevant_solution;
I proceed by reading two separate external meshes and attaching them to the
appropriate triangulations.
Then I begin the Setup_System(). I start by getting the number of DOFs
from each of the triangulations and adding them
in a vector of the DOFs for each block. I am saying that the Solid is
block 0 and the Shell is block 1.
const unsigned int n_dof_sld = dof_handler_sld.n_dofs(); // = 24
const unsigned int n_dof_sh = dof_handler_sh.n_dofs();// = 48
const std::vector dofs_per_block = {n_dof_sld, n_dof_sh};
I then begin work on the sparsity pattern:
BlockDynamicSparsityPattern dsp(dofs_per_block , dofs_per_block);
for (unsigned int i=0; iFrom what I can see, the system_matrix looks right with n_blocks = 2 and
start_indices = {0, 24, 72}, as expected.
Next, I run through the Assemble_System() to create a local cell matrix. I
start with the Solid elements and cell_matrix_sld. I'm skipping the
calculation of the cell matrix for brevity's sake.
for (const auto : dof_handler_sld.active_cell_iterators())
{
// *** Placeholder for calculating cell_matrix_sld, a 24x24 matrix ***
cell->get_dof_indices(local_dof_indices_sld);
// Add the cell matrix to the system matrix
for (unsigned int i = 0; i < dofs_per_cell_sld; ++i)
for (unsigned int j = 0; j < dofs_per_cell_sld; ++j)
system_matrix.add(local_dof_indices_sld[i],
local_dof_indices_sld[j],
cell_matrix_sld(i, j));
} // end of cell loop
The bit of code 'system_matrix.add(...)' works for when i=0 and j=0, but
fails on the next loop when j=1. I get this error:
An error occurred in line <1931> of file <.../deal.II/lac/sparse_matrix.h>
in function
void
dealii::SparseMatrix::add(dealii::SparseMatrix::size_type,
dealii::SparseMatrix::size_type, number) [with number = double;
dealii::SparseMatrix::size_type = unsigned int]
The violated condition was:
(index != SparsityPattern::invalid_entry) || (value == number())
Additional information:
You are trying to access the matrix entry with index <0,1>, but this
entry does not exist in the sparsity pattern of this matrix.
The most common cause for this problem is that you used a method to
build the sparsity pattern that did not (completely) take into account
all of the entries you will later try to write into. An example would
be building a sparsity pattern that does not include the entries you
will write into due to constraints on degrees of freedom such as
hanging nodes or periodic boundary conditions. In such cases, building
the sparsity pattern will succeed, but you will get errors such as the
current one at one point or other when trying to write into the
entries of the matrix.
So, it appears that I'm not setting something up correctly with the
system_matrix and sparsity pattern. Step-20 uses the same
system_matrix.add(...) approach, and does not have any additional
preparation of the system_matrix that I can see.
One thing I skipped is this line from Step 20:
DoFTools::make_sparsity_pattern
(dof_handler,
dsp);
I'm unable to use this line because I have the two separate DOF handlers.
Do you have any idea of what I can do here?
Best regards,
Alex
On Saturday, October 28, 2023 at 6:58:21 PM UTC-4 Wolfgang Bangerth wrote:
> On 10/27/23 12:21, Daniel Arndt wrote:
> >
> > If you can
Re: [deal.II] Mixed codimensionality
Hi Daniel, Thanks for the response. I am essentially using the approach in Step-8, where I am seeking a vector-valued displacement solution using FE_System. So, I've been using: dealii::FESystem<3>(FE_Q<3>(1), 3 ) and dealii::FESystem<2,3>(FE_Q<2,3>(1), 3 ) Based on your suggestion, I am picturing doing something like this: // Shared Objects AffineConstraints constraints; SparsityPatternsparsity_pattern; SparseMatrix system_matrix; Vector solution;// This is the solution vector Vector system_rhs; // This is the load / right-hand-side vector // Dimension dependent objects // 3D Triangulation<3,3> triangulation_3D; DoFHandler<3,3> dof_handler_3D; FESystem<3,3>fe_3D; // 2D Triangulation<2,3> triangulation_2D; DoFHandler<2,3> dof_handler_2D; FESystem<2,3>fe_2D; // 1D Triangulation<1,3> triangulation_1D; DoFHandler<1,3> dof_handler_1D; FESystem<1,3>fe_1D; Then I need to proceed with each of these triangulations, and *somehow* link them together within the system_matrix and system_rhs. I have been building the system matrix using this line for each cell in the dof_handler.active_cell_iterators: constraints.distribute_local_to_global( cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); I see two challenges using this: 1) keeping track of the dof_indices across 3 dof_handlers so they can be added to the correct place in the system_matrix 2) Connecting the triangulations within the system matrix. Have I understood your suggestion correctly? Do you forsee any other challenges? Thanks very much, Alex -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/7b15f86c-0069-42c9-b7ed-4bde752ee07bn%40googlegroups.com.
[deal.II] Mixed codimensionality
Dear all, I've been working with deal.ii in the solid mechanics field for both 3D and 2D cases. I am now being asked by my colleagues to develop a model that can combine elements of different dimensionality. This is approach is already implemented in my colleagues' abaqus models, and so we're looking to replicate the functionality in deal.ii. The elements to be combined would be FE_System elements of type: - 3D volumetric element in 3D space - 2D planar element in 3D space - 1D linear element in 3D space I've attached a sketch of an application that might use this. >From what I've seen, deal.ii does not like to mix 's. Does anyone have thoughts on how I might go about implementing a model with these elements? If I need to modify or implement new functionalities, could you estimate the magnitude of time required to make changes? (i.e. days, weeks, months) I realize that there are alternative methods to represent the physical system that could avoid this problem. However, I am essentially being asked to automate the conversion from existing abaqus input files into deal.ii. I appreciate any ideas you all might have. Thanks, Alex -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/a7ff55e9-b406-403b-b176-75690c2d121en%40googlegroups.com.
Re: [deal.II] Non-polynomial Shape Functions
Thanks, Wolfgang. Your solution worked for me.
Just for posterity, I got stuck again because I created FE_System and
FE_Enriched in the wrong order:
const FESystem fesys;
const FE_Enriched fenr;
my problem was resolved when I switched it around to:
const FE_Enriched fenr;
const FESystem fesys;
Cheers,
Alex
On Wednesday, June 21, 2023 at 1:25:44 AM UTC-4 Wolfgang Bangerth wrote:
> On 6/20/23 10:00, Alex Quinlan wrote:
> >
> > How can I go about getting a vector solution with FE_Enriched? I have
> tried
> > to use vector arguments in FE_Enriched for base elements, enriched
> elements,
> > and enrichment functions, based on the documentation. I have also tried
> to
> > use FE_Enriched as an argument to FE_System:
> >
> > const FESystem fesys;
> > const FE_Enriched fenr;
> >
> > shElast::shElast()
> > :fesys(fenr(FE_Q(2),
> > FE_Q(1),
> >)
> > , 3)
>
> This is the right idea, but the wrong syntax. In constructors, what comes
> after the : is a comma-separated list of the form
> member_var_name(arg1,arg2,...),
> ...
>
> So you need to write this list as
>
> shElast::shElast()
> :fenr(FE_Q(2),
> FE_Q(1),
> )
> ,fesys(fenr, 3)
> {...}
>
> Best
> W.
> --
>
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/d00a28c2-5053-4545-8aa6-0059fa274bd5n%40googlegroups.com.
Re: [deal.II] Non-polynomial Shape Functions
Hi Wolfgang, Thanks for the guidance on the pointer. I have gotten a compiled program using: const FE_Enriched fenr; PotentialFunction potential; shElast::shElast() :fenr( FE_Q(2), FE_Q(1), ) However, I am in need of a 3D vector solution. Previously, I have successfully used FE_System for non-enriched 3D vector problems, i.e.: const FESystem fesys; shElast::shElast() :fesys( FE_Q(2), 3) How can I go about getting a vector solution with FE_Enriched? I have tried to use vector arguments in FE_Enriched for base elements, enriched elements, and enrichment functions, based on the documentation. I have also tried to use FE_Enriched as an argument to FE_System: const FESystem fesys; const FE_Enriched fenr; shElast::shElast() :fesys( fenr( FE_Q(2), FE_Q(1), ) , 3) Are either of these directions valid? I can't tell if I'm failing to use the correct syntax or if I'm completely off base with the approach. Thanks, Alex -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/1e642982-e834-47e0-b608-360d12a96bd3n%40googlegroups.com.
Re: [deal.II] Non-polynomial Shape Functions
Hi Wolfgang,
I've been looking at a few sources regarding FE_Enriched:
(A) https://www.dealii.org/current/doxygen/deal.II/classFE__Enriched.html
(B)
https://github.com/dealii/dealii/blob/master/tests/fe/fe_enriched_step-36b.cc
I've been trying to get my problem class to work with a placeholder
enrichment function from example (B), but I'm getting tripped up.
I looked at the formatting on Line 146 of fe_enriched.cc:
FE_Enriched::FE_Enriched(
const FiniteElement _base,
const FiniteElement _enriched,
const Function * enrichment_function)
I have tried to implement this format in the constructor for my problem
class. I am trying to set fe_base as FE_Q, and then fe_enriched as another
FE_Q multiplied
by the enrichment function.
template
shElast::shElast(
const std::string & initial_mesh_filename,
const std::string & output_filename)
:
*fe( FE_Q(1) , FE_Q(1) , enrichment)*
, dof_handler(triangulation)
, quadrature_formula(fe.degree + 1)
, initial_mesh_filename(initial_mesh_filename)
, output_filename(output_filename)
, mapping( FE_Q(1))
{}
where I've previously declared:
Triangulation triangulation;
* const FE_Enriched fe;*
DoFHandlerdof_handler;
const QGauss quadrature_formula;
const MappingFE mapping;
* EnrichmentFunctionenrichment;*
*When I try to compile my program, I get this error:*
dealii/examples/enrich/enrich9.cc:324:25: error: no matching function for
call to ‘dealii::FE_Enriched<3, 3>::FE_Enriched(dealii::FE_Q<3, 3>,
dealii::FE_Q<3, 3>, problem_namespace::EnrichmentFunction<3>&)’
324 | , mapping( FE_Q(1))
*Then I get the following notes about my arguments and failure to find
conversions:*
In file included from dealii/examples/enrich/enrich9.cc:53:
fe/fe_enriched.h:263:3: note: candidate: ‘dealii::FE_Enriched::FE_Enriched(const std::vector*>&, const
std::vector&, const
std::vector*(const typename dealii::Triangulation::cell_iterator&)> > >&) [with int dim = 3; int spacedim = 3;
typename dealii::Triangulation::cell_iterator =
dealii::TriaIterator >]’
263 | FE_Enriched(
| ^~~
deal.II/fe/fe_enriched.h:264:62: note: no known conversion for argument 1
from ‘dealii::FE_Q<3, 3>’ to ‘const std::vector*, std::allocator*> >&’
264 | const std::vector *> ,
| ~^~~
deal.II/fe/fe_enriched.h:249:3: note: candidate: ‘dealii::FE_Enriched::FE_Enriched(const dealii::FiniteElement*, const std::vector*>&, const
std::vector*(const typename dealii::Triangulation::cell_iterator&)> > >&) [with int dim = 3; int spacedim = 3;
typename dealii::Triangulation::cell_iterator =
dealii::TriaIterator >]’
249 | FE_Enriched(
| ^~~
deal.II/fe/fe_enriched.h:250:62: note: no known conversion for argument 1
from ‘dealii::FE_Q<3, 3>’ to ‘const dealii::FiniteElement<3, 3>*’
250 | const FiniteElement *
fe_base,
| ~^~~
deal.II/fe/fe_enriched.h:200:3: note: candidate: ‘dealii::FE_Enriched::FE_Enriched(const dealii::FiniteElement&) [with int dim = 3; int spacedim = 3]’
200 | FE_Enriched(const FiniteElement _base);
| ^~~
deal.II/fe/fe_enriched.h:200:3: note: candidate expects 1 argument, 3
provided
deal.II/fe/fe_enriched.h:186:3: note: candidate: ‘dealii::FE_Enriched::FE_Enriched(const dealii::FiniteElement&, const dealii::FiniteElement&, const dealii::Function*) [with int dim = 3;
int spacedim = 3]’
186 | FE_Enriched(const FiniteElement _base,
| ^~~
deal.II/fe/fe_enriched.h:188:51: note: no known conversion for argument 3
from ‘shell::EnrichmentFunction<3>’ to ‘const dealii::Function<3, double>*’
188 | const Function *
enrichment_function);
|
^~~
It seems that I'm misunderstanding a few things about how FE_Enriched
works. Can you see what I'm doing wrong in constructing the FE_Enriched?
It looks like there's issues with forms I am using for FE_Q and the
enrichment function. Do you have any tips or suggested resources?
Many thanks,
Alex
On Friday, June 9, 2023 at 8:27:37 AM UTC-4 Alex Quinlan wrote:
> Thanks, Wolfgang. I'll dig into FE_Enriched!
>
> On Thursday, June 8, 2023 at 6:57:44 PM UTC-4 Wolfgang Bangerth wrote:
>
>> On 6/8/23 09:40, Alex Quinlan wrote:
>> >
>> > I am looking to modify the FE_Q element to have non-polynomial shape
>> > functions. It looks like FE_Q_Base takes an argument of type
>> > ScalarPolynomialBase (fe_q_base.cc: line 418).
>> >
>> > I am looking at a 2D ele
Re: [deal.II] Non-polynomial Shape Functions
Thanks, Wolfgang. I'll dig into FE_Enriched! On Thursday, June 8, 2023 at 6:57:44 PM UTC-4 Wolfgang Bangerth wrote: > On 6/8/23 09:40, Alex Quinlan wrote: > > > > I am looking to modify the FE_Q element to have non-polynomial shape > > functions. It looks like FE_Q_Base takes an argument of type > > ScalarPolynomialBase (fe_q_base.cc: line 418). > > > > I am looking at a 2D element and instead of 2D polynomial shape > functions like > > this: > > > > N = (1 - ξ)(1 - η) > > > > I would like to add some additional terms like this: > > > > N = (1 - ξ)(1 - η) + ζ > > > > where ζ changes at different quadrature points. > > > > It seems that I would need to access the values of ζ before computing > > fe_values, since the shape functions and their derivatives would be > dependent > > on ζ. > > > > Any thoughts on how feasible this would be to implement? > > Alex: > if your shape functions are not polynomial, then this is the wrong > approach -- > FE_Q_Base is derived from FE_Poly, which is specifically built to create > polynomial bases. > > I'm not sure we have a particularly good starting point for a completely > general choice of ζ(ξ,η), but you might want to take a look at what > FE_Enriched is doing. > > Best > W. > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/97405fd5-9732-4fb2-b1b7-ce6f2d1a155en%40googlegroups.com.
[deal.II] Non-polynomial Shape Functions
Hi everyone, I am looking to modify the FE_Q element to have non-polynomial shape functions. It looks like FE_Q_Base takes an argument of type ScalarPolynomialBase (fe_q_base.cc: line 418). I am looking at a 2D element and instead of 2D polynomial shape functions like this: N = (1 - ξ)(1 - η) I would like to add some additional terms like this: N = (1 - ξ)(1 - η) + ζ where ζ changes at different quadrature points. It seems that I would need to access the values of ζ before computing fe_values, since the shape functions and their derivatives would be dependent on ζ. Any thoughts on how feasible this would be to implement? Thanks for your input, Alex -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/0671820a-700c-4bbc-98f3-ba53ea5c7b76n%40googlegroups.com.
Re: [deal.II] Importing nodal BCs and accessing vertices
Thanks Wolfgang. I appreciate the feedback. I'll see if I can implement some of the cost-saving tips that you suggested. -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/2c35e77b-3c66-4683-8186-6951b103142cn%40googlegroups.com.
[deal.II] Importing nodal BCs and accessing vertices
Hi everyone,
I am working on adapting some abaqus jobs for deal.ii. In the abaqus
format, It is typical to specify boundary conditions in terms of (node
number, dof number, magnitude).
I have found a working approach for doing this, but I am looking to find a
more efficient way as I may end up with millions of nodes and thousands of
vertices. Here's an example of how I have been applying nodal displacement
BCs:
// define a "node" by the coordinates
const Point proot1(0., 0.0209, -0.0101);
// loop through all cells
for (const auto : dof_handler.active_cell_iterators())
{
// loop through all of the vertices in a given cell
for (unsigned int v = 0; vn_vertices() ; v++)
{
// check if a vertex is near to the node.
if (cell->vertex(v).distance(proot1) < 0.001) //
{
// apply a constraint
constraints.add_line(cell->vertex_dof_index(v,0, cell->active_fe_index()
));
} } }
In an ideal world, I could compute the global dof number based on the node
number and local dof number. Something like:
global_dof = node_number * 3 + local_dof
constraints.add_line(global_dof)
However, I understand that this approach isn't implemented in deal.ii
because it would be useless after mesh refinement.
Is there a good way to find the global_dof numbers for a vertex/node?
If not, can I iterate through all of the vertices directly (rather than
going through the cell)? Something like:
for (const auto : dof_handler.vertex_iterators())
{
if (vertex.distance(proot1) < 0.001)
{
constraints.add_line(vertex->vertex_dof_index(0) );
}}
I'm interested on your thoughts on how to efficiently import nodal boundary
conditions from an external file. I couldn't find any other posts on this
specific topic, but I apologize if this is a duplicate.
Thanks,
Alex
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/7a0ff20d-7d25-4b2d-a094-5bccc1388a5fn%40googlegroups.com.
Re: [deal.II] 3D Mixed Simplex Mesh
Dear Wolfgang,
Thanks very much for your reply. I'm happy to hear that you think my
problems can be resolved without too much difficulty. I'd like to try
writing a patch for the mixed 3D meshes. I have done some development work
before, but I have never sent it out of our group; your advice on writing
the steps to write a patch would be most welcome!
I will re-review the relevant code and see if I can put together a general
plan on how to approach this.
Best regards,
Alex
On Wednesday, November 9, 2022 at 4:44:25 PM UTC-5 Wolfgang Bangerth wrote:
> Alex:
>
> > I am coming from a background of solid mechanics FEA and I've been using
> > CalculiX as an openSource program for some time. My group has decided
> > to switch to deal.ii
>
> Good choice -- we applaud you :-)
>
>
> > Part 1: Abaqus .inp format
> > I had success importing hex meshes from the abaqus <.inp> format, but it
> > seems that the UCD reader does not accept tet meshes. I am basing this
> > on the following bit of code in grid_in::read_ucd around Line 948
> >
> > if (((cell_type == "line") && (dim == 1)) ||
> > ((cell_type == "quad") && (dim == 2)) ||
> > ((cell_type == "hex") && (dim == 3)))
> > // found a cell
> > {
> >
> > The absence of
> >((cell_type == "tet") && (dim == 3))
> > suggests that simplex elements are not supported when importing from
> > .inp. Am I correct in this conclusion. How difficult would it be for
> > me to add this functionality in?
>
> Not very difficult in all likelihood. The key difficulty in adding
> support for other cells in this function is that it uses the
> GeometryInfo class, which only provides information about hypercube
> cells. For example, we read vertex indices like this in that function:
> for (const unsigned int i : GeometryInfo::vertex_indices())
> in >> cells.back().vertices[GeometryInfo::ucd_to_deal[i]];
>
> This needs to be changed to account for the number of vertices a cell
> has, given the cell type, not the dimension. The right approach will be
> to use the ReferenceCell objects defined in namespace ReferenceCells
> (plural), and the other readers in that function will serve as guidance.
> You might want to take a look at the gmsh reader, for example.
>
>
> > Part 2: GMSH format
> > While looking at grid_in.cc, I saw that the MSH reader supported tet
> > elements. I found the FAQ page supplied these examples of simplex
> > meshes:
> https://www.dealii.org/developer/doxygen/deal.II/group__simplex.html
> >
> > I first ran the 2D simplex mesh example, and then I modified it to work
> > for 3D without much problem. Then, I switched to the mixed mesh, and
> > had success running the 2D example. However, I started having problems
> > when I switched to 3D. The mesh seems to import, but then I have issues
> > operating with it. I found several places within fe_simplex_p.cc (L714,
> > L759) where the dimension is asserted to be 2:
> >
> >AssertDimension(dim, 2);
> > This seems like a problem for a 3D analysis, so I am wondering if I need
> > to do something different so that this function doesn't get called.
>
> No, I think this is just not implemented. The function is called on
> mixed meshes when the library wants to know how elements of type
> FE_SimplexP on tetrahedra connect to another element types on other cell
> shapes (namely, on wedges and pyramids). The comparison is simply not
> implemented, but if you want to compute on mixed 3d meshes, you will
> have to implement the missing parts of this function. As long as you
> want to use the same polynomial degree on all cells, this should not
> even be particularly difficult and I would be happy to walk you through
> the necessary steps if you were interested in writing up a patch!
>
> Best
> Wolfgang
>
> --
>
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/885a1c8c-a301-4467-8d16-9a5d156d6f95n%40googlegroups.com.
[deal.II] Re: 3D Mixed Simplex Mesh
I realized that the mesh file I provided had an issue. The corner nodes of
the tetrahedral elements coincided with a mid-face node of the hexahedron
element. I've fixed that and uploaded a new mesh.
Thanks,
Alex
On Monday, November 7, 2022 at 3:25:25 PM UTC-5 Alex Quinlan wrote:
> Hello everyone,
> I have been reading the messages in this user group for a while, but this
> is my first post! I am hoping that you can help with my problem.
>
> I am coming from a background of solid mechanics FEA and I've been using
> CalculiX as an openSource program for some time. My group has decided to
> switch to deal.ii to implement some advanced functionality that the other
> programs don't have. I'm an engineer with some programming experience,
> though fairly new to C++.
>
> We are working with some complex geometries that are not practical to mesh
> with only hex elements. The meshes are generated externally, so I am
> trying to import the meshes into deal.ii. I have tried to do this in two
> general ways and run into problems with each, so I will split my post into
> 2 parts.
>
> Part 1: Abaqus .inp format
> I had success importing hex meshes from the abaqus <.inp> format, but it
> seems that the UCD reader does not accept tet meshes. I am basing this on
> the following bit of code in grid_in::read_ucd around Line 948
>
> if (((cell_type == "line") && (dim == 1)) ||
> ((cell_type == "quad") && (dim == 2)) ||
> ((cell_type == "hex") && (dim == 3)))
> // found a cell
> {
>
> The absence of
>((cell_type == "tet") && (dim == 3))
> suggests that simplex elements are not supported when importing from
> .inp. Am I correct in this conclusion. How difficult would it be for me
> to add this functionality in?
>
>
> Part 2: GMSH format
> While looking at grid_in.cc, I saw that the MSH reader supported tet
> elements. I found the FAQ page supplied these examples of simplex meshes:
> https://www.dealii.org/developer/doxygen/deal.II/group__simplex.html
>
> I first ran the 2D simplex mesh example, and then I modified it to work
> for 3D without much problem. Then, I switched to the mixed mesh, and had
> success running the 2D example. However, I started having problems when I
> switched to 3D. The mesh seems to import, but then I have issues operating
> with it. I found several places within fe_simplex_p.cc (L714, L759) where
> the dimension is asserted to be 2:
>
>AssertDimension(dim, 2);
>
> This seems like a problem for a 3D analysis, so I am wondering if I need
> to do something different so that this function doesn't get called.
>
> The problem seems to arrive when I call dof_handler.distribute_dofs(fe);
> I've attached my modified version of the example program along with my
> simple mixed-element 3D mesh.
>
> Is there something that I am missing regarding the dof_handler when trying
> to use a 3D mixed-element mesh?
>
>
> Thanks very much,
> Alex
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/83b8a49e-17e1-4f5b-aafa-4a6438bb1073n%40googlegroups.com.
mixed2top.msh
Description: Mesh model
[deal.II] 3D Mixed Simplex Mesh
Hello everyone,
I have been reading the messages in this user group for a while, but this
is my first post! I am hoping that you can help with my problem.
I am coming from a background of solid mechanics FEA and I've been using
CalculiX as an openSource program for some time. My group has decided to
switch to deal.ii to implement some advanced functionality that the other
programs don't have. I'm an engineer with some programming experience,
though fairly new to C++.
We are working with some complex geometries that are not practical to mesh
with only hex elements. The meshes are generated externally, so I am
trying to import the meshes into deal.ii. I have tried to do this in two
general ways and run into problems with each, so I will split my post into
2 parts.
Part 1: Abaqus .inp format
I had success importing hex meshes from the abaqus <.inp> format, but it
seems that the UCD reader does not accept tet meshes. I am basing this on
the following bit of code in grid_in::read_ucd around Line 948
if (((cell_type == "line") && (dim == 1)) ||
((cell_type == "quad") && (dim == 2)) ||
((cell_type == "hex") && (dim == 3)))
// found a cell
{
The absence of
((cell_type == "tet") && (dim == 3))
suggests that simplex elements are not supported when importing from .inp.
Am I correct in this conclusion. How difficult would it be for me to add
this functionality in?
Part 2: GMSH format
While looking at grid_in.cc, I saw that the MSH reader supported tet
elements. I found the FAQ page supplied these examples of simplex meshes:
https://www.dealii.org/developer/doxygen/deal.II/group__simplex.html
I first ran the 2D simplex mesh example, and then I modified it to work for
3D without much problem. Then, I switched to the mixed mesh, and had
success running the 2D example. However, I started having problems when I
switched to 3D. The mesh seems to import, but then I have issues operating
with it. I found several places within fe_simplex_p.cc (L714, L759) where
the dimension is asserted to be 2:
AssertDimension(dim, 2);
This seems like a problem for a 3D analysis, so I am wondering if I need to
do something different so that this function doesn't get called.
The problem seems to arrive when I call dof_handler.distribute_dofs(fe);
I've attached my modified version of the example program along with my
simple mixed-element 3D mesh.
Is there something that I am missing regarding the dof_handler when trying
to use a 3D mixed-element mesh?
Thanks very much,
Alex
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/a905250a-c35f-4a75-952a-a2237451e94fn%40googlegroups.com.
/* -
*
* Copyright (C) 1999 - 2021 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* -
*/
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
using namespace dealii;
class Step3
{
public:
Step3();
void run();
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
Triangulation<3> triangulation;
const hp::MappingCollection<3> mapping;
const hp::FECollection<3> fe;
const hp::QCollection<3> quadrature_formula;
DoFHandler<3> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix system_matrix;
Vector solution;
Vector system_rhs;
};
Step3::Step3()
: mapping(MappingFE<3>(FE_SimplexP<3>(1)), MappingFE<3>(FE_Q<3>(1)))
, fe(FE_SimplexP<3>(2), FE_Q<3>(2))
, quadrature_formula(QGaussSimplex<3>(3), QGauss<3>(3))
, dof_handler(triangulation)
{}
void Step3::make_grid()
{
GridIn<3>(triangulation).read("input/mixed.msh");
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
void Step3::setup_system()
{
