Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Linear constraints #401

Merged
merged 1 commit into from
Feb 7, 2022
Merged

Linear constraints #401

merged 1 commit into from
Feb 7, 2022

Conversation

lijas
Copy link
Collaborator

@lijas lijas commented Nov 21, 2021

This wip PR makes it possible to use linear constraints. I was able to use a lot of the old structure of the ConstraintHandler which is nice. I also got inspiration on how to condense the stiffness matrix from deal ii :)

Currently this PR is not very fail safe, meaning that the user can add Dirichlet and/or linear constraints that does not "work" together, without warning the user. There is no logic that assures that constraints are linearly independent etc (this is something the user must make sure themselves).

When adding constraints, the user must reference global dof ids directly. For this, we need better tools for querying dof ids on faces/edges/nodes etc.

Example:

...
ch = ConstraintHandler(dh)
add!(ch, Dirichlet(:u, getfaceset(grid, "left"), (x,t)->0.0))
    add!(ch, LinearConstraint(4, [(7 => 2.0)], 1.5)) # u4 = 2.0*u7 - 1.5
close!(ch)
update!(ch, 0.0)
...
apply!(K, f, ch)
a = K\f
apply!(a,ch)

@lijas lijas mentioned this pull request Nov 26, 2021
@termi-official
Copy link
Member

Looks cool! Just a quick question here again, anything wrong with using the maps returned by __close! to identify the dofs?

@lijas
Copy link
Collaborator Author

lijas commented Nov 29, 2021

For node dofs that should work fine i think.

But regarding the dofs, I am thinking about adding some dof-tools to get dof-ids on faces/nodes.

@kimauth
Copy link
Member

kimauth commented Nov 30, 2021

Maybe some refactoring of the code in the ConstraintHandler can help for the dof-tools? At least the ConstraintHandler is able to find the dofs related to faces and also to certain node numbers.

@koehlerson
Copy link
Member

given a face and/or edge you should still be able to figure out the corresponding dofs with those dicts https://github.com/Ferrite-FEM/Ferrite.jl/blob/master/src/Dofs/DofHandler.jl#L281 But we throw them away everytime we call close! instead of __close!

@lijas
Copy link
Collaborator Author

lijas commented Nov 30, 2021

I dont think the facedict from close! can be used... It stores (thee) nodeids of the faces, but usually we references faces via FaceIndex (cellid+localfaceid).

I think that #360 would be a better approach (and relatively efficient). (facedofs!(dofvec, dh, FaceIndex(cellid,lid)))

@koehlerson
Copy link
Member

yes, but with sortface(face) you get the tuple you need (at least I think so). FaceIndex recomputation to global sorted face is also always possible via sorface(faces(grid.cells[faceidx[1])faceidx[2]))

@codecov-commenter
Copy link

codecov-commenter commented Nov 30, 2021

Codecov Report

Merging #401 (cc88da3) into master (636f273) will increase coverage by 1.48%.
The diff coverage is 95.18%.

❗ Current head cc88da3 differs from pull request most recent head b82cda1. Consider uploading reports for the commit b82cda1 to get more accurate results

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #401      +/-   ##
==========================================
+ Coverage   89.78%   91.27%   +1.48%     
==========================================
  Files          22       22              
  Lines        2701     3036     +335     
==========================================
+ Hits         2425     2771     +346     
+ Misses        276      265      -11     
Impacted Files Coverage Δ
src/Dofs/MixedDofHandler.jl 86.36% <66.66%> (-0.61%) ⬇️
src/Dofs/DofHandler.jl 83.39% <77.77%> (-2.71%) ⬇️
src/Dofs/ConstraintHandler.jl 94.62% <97.09%> (+1.17%) ⬆️
src/interpolations.jl 86.84% <0.00%> (-3.99%) ⬇️
src/Export/VTK.jl 91.22% <0.00%> (-1.20%) ⬇️
src/FEValues/common_values.jl 88.72% <0.00%> (-0.17%) ⬇️
src/Grid/grid.jl 86.62% <0.00%> (ø)
src/L2_projection.jl 100.00% <0.00%> (ø)
src/FEValues/cell_values.jl 100.00% <0.00%> (ø)
... and 4 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 636f273...b82cda1. Read the comment docs.

@lijas
Copy link
Collaborator Author

lijas commented Dec 2, 2021

I found a bug where I was pushing new entries to K while thinking I was only reading the old entries

for r in nzrange(K, col)
   row = K.rowval[r]
   val = K.nzval[r]
   ...
   K[row, col] += val #Can possibly create a new entry i K which modifies K.rowval[]

so I had to change it to

range = nzrange(K, col)
_rows = K.rowval[range]
_vals = K.nzval[range]
for (row, val) in zip(_rows, _vals)
   ...
   K[row, col] += val

but this allocates _rows and _vals which is not super nice.

We could force the user to first call create_sparsity_pattern(dh, ch), which means that we could go with the first version (because K[row, val] += val would not add any new entries in K). However, there is no way to check that K actually has the correct sparsity-pattern needed for the condensation process (i.e they might have called create_sparsity_pattern(dh) without ch). This means that the user might call apply!(K,f, ch) with the wrong sparsity-pattern, and not get any error message...

Thoughts on what we should do?

Edit: the assembler class also requires the user to have called create_sparsity_pattern(dh), but at-least it throws an error if they have not done it.

@kimauth
Copy link
Member

kimauth commented Dec 2, 2021

Can't we check if K has non-zero (well or explicit zero) values in these positions and at least throw a warning when a user hasn't created the right sparsity pattern?

@termi-official
Copy link
Member

I think the simplest approach here is to write a new function create_condensed_sparsity_pattern(dh,ch) or maybe a dispatch create_sparsity_pattern(dh,ch) which creates the correct pattern for the condensed matrix. If you want to work with direct solvers or if the matrix needs to be represented explicitly.

If just iterative solvers are used, then maybe just computing the action of the "condensation matrix" C (and C') could be more efficient - at least in some cases.

@lijas
Copy link
Collaborator Author

lijas commented Dec 3, 2021

@kimauth I think that is pretty much the same amount of work as doing the condensation.
But this gave me an idea, the K[row, col] += val operation probably do some internal checking if [row, col] entry exists, so maybe I can use that in some way.

@termi-official That is what I have added already :) create_sparsity_pattern(dh,ch). The thing is, if they do not call that function (for example only calling create_sparsity_pattern(dh)), and then try to condense the system, they first code snippet I posted will condense the system wrongly, and will not throw an error.

@termi-official
Copy link
Member

Ah, sorry, no idea how I missed that in the discussion.

@lijas
Copy link
Collaborator Author

lijas commented Jan 25, 2022

I have tested this in my own code and everything seems to work.
I compared the performance of using linear constraints to normal Dirichlet constraints, and it seems like they have comparable times when condensing the stiffness matrix (linear constraints will not be a major bottleneck in the code at-least).

My next TODO is to add some docs, but if anyone wants to start reviewing this PR, please feel free too :P (@fredrikekre)

Copy link
Member

@termi-official termi-official left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks really cool! I left some comments on the docs.

Ping @koehlerson we could utilize this for AMR, as this is basically just a different approach to applying the P-matrix, which we implemented for the AMR student project.

docs/src/manual/constraints.md Outdated Show resolved Hide resolved
docs/src/manual/constraints.md Outdated Show resolved Hide resolved
@fredrikekre
Copy link
Member

Very cool.

@lijas
Copy link
Collaborator Author

lijas commented Feb 7, 2022

Nice, thanks for review and updates.
I will update the docs in another PR

koehlerson pushed a commit that referenced this pull request Apr 22, 2022
This adds the AffineConstraint constraint type that can be added to the
constraint handler. An affine constraint of the form
    u_i = \sum_j (c_j * u_j) + b,
where u_i are dofs, c_j scaling factors and b a inhomogeneity is
constructed and added to the constraint handler as
    ac = AffineConstraint(i, [u_j => c_j, ...], b)
    add!(ch, ac)

Sparsity patterns must then be constructed with the constraint handler
e.g. create_sparsity_pattern(dh, ch), and apply! must be used after
solving the system (apply!(u, ch), where u is the solution).
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

6 participants