| Thread Previous • Date Previous • Date Next • Thread Next |
Johan Hoffman wrote:
Patrick Riesen wrote:We apply this type of boundary conditions to the discrete matrix equation after having assembled the system, see the implementation in Unicorn. This is not very neat, and I think you instead probably would like to apply strong bc through the definition of the function space.Johan Hoffman wrote:Patrick Riesen wrote:Patrick Riesen wrote:Johan Hoffman wrote:Patrick Riesen wrote:If your boundary is curved (that is: not a line or plane), applying your normal bc (slip bc) weakly in this way based on the default discontinuous face/edge normals will result in artificial friction, or if your penalty is large enough zero velocity.Patrick Riesen wrote:Garth N. Wells wrote:On 17/05/10 08:50, Patrick Riesen wrote:Anders Logg wrote:On Fri, May 07, 2010 at 02:39:41PM +0200, Patrick Riesen wrote:hello, i want to zero the normal velocity on a boundary part of my 2D-domain, so i used something like n = cell.n a = ... + inner(w,dot(v,n)*n)*ds(2) + ... in my form-file (w=testfunction, v=trialfunction). i marked the boundary part using the "exterior facet domains" meshfunction, which will be picked up in assemble(..): MeshFunction< uint >* _ext_facet_domains = mesh.data().create_mesh_function("exterior facet domains", 1) *_ext_facet_domains = 5; boundary.mark(*_ext_facet_domains, 2)can i do this? or what is wrong with the above way? because i don'tget zero velocities... thank your for help, patrickWhat if you add a (big) parameter in front of that term? Does it haveany effect at all on the solution? I would also suggest writing the term as dot(w, n)*dot(v, n)*ds(2)thank you for the help anders, so i noticed the following:without parameter in front, i could see no effect at all. if i add a parameter C in front and increase it stepwise until C is about 10^6 ican observe how the velocities cancel out on that boundary.You can probably improve this significantly by enforcing zero normal velocity using Nitsches's method. This involves just adding two extra terms to your form. See demo/pde/dg/poisson for an example.Garthhello, so based on the literature from peter hansbo http://www.math.chalmers.se/~hansbo/nitscherev_preprint.pdfi used the nitsche's method so impose a ALE boundary condition of the form(v - v_m)*n = 0 on my (moving) boundary (garth: it ends up to 5 additional terms). (v=fluid velocity, v_m = mesh velocity, n = outward normal).I think it works correctly but i observe a loss of volume because the above bcs is only fullfilled to differences to about 1.0e-4 to 1.0e-5, so in every timestep i loose some mass/volume due to that, i guess.in hansbo's script there is not much about the stability parameter h^-1 needed and i just used the same as in the dg/poisson example, but does somebody know more about how to choose or estimate it?regards, patrickyes, i noticed when my penalty coefficient is large, it doesn't help the solution, but i just get closer to v = v_m.What you can do is:(1) apply your bc strongly at each surface node instead, based on a vertex normal, e.g. from a weighted average from surrounding face/edge normals (this is also mass conservative), see:this was where i first started before the weak imposition, but i had difficulties in creating those vertex normals.Ok. You can check in the Unicorn solver, this is what we use: http://www.fenics.org/wiki/Unicornhttp://www.csc.kth.se/~jhoffman/archive/papers/scm.pdf(2) form a continuous normal field to apply your bc weakly, for example a p.w. linear normal field over the surface. This may also work depending on what method you use to weakly apply your bc.so i would solve for a normal field? can you elaborate some more on this?do you think i can keep the weak statement and just replace the facet normal 'n' (the facet normal from ufl) by a linear normal function of the vertex normals?Yes, I would expect this to work.ok, i think it did not work (or not much better). I might try the strong way.But how to specify strong Dirichlet conditions for a variable 'v_m' as a function of an other variable 'v' in the mixed problem (v_m, v) in fenics?
so far i came up with something similar. i use the dirichletBC class together with a special expression which sets the vertex values on the free surface to v_m = v in the space of v_m, i.e the boundary is fully lagrangian. This is more restrictive than specifying the condition on the normal components, but was easier to implement for initial testing.
i here plotted the element-wise area differences across timestepping: http://people.ee.ethz.ch/~priesen/public/tmsoim_ale.gifi still get a total cumulative area loss of about 10e-2 from the O(1) square. I'm not sure if this is acceptable. although, the deformation is
rather significant. should i be happy or not? ;-)an other approach on the formulation level which might be possible is to use the Lifting functionality in ufl for the dirichlet moving boundary part. i will look into that.
thanks for the help so far, regards, patrick
Best, Johanregards, patrickBest, JohanGood luck, Johanthank you (for the help) ;-) regards, patricksounds interesting, i'll look into it. thank you, patricki guess this approach works, but there is something else:looking at the (time-dependent) solution (*pvd file set or single vtufiles) in paraview, i ran into several ERROR: In/home/kitware/Dashboard/MyTests/ParaView-3-8/ParaView-3.8/ParaView/VTK/IO/vtkXMLDataReader.cxx,line 509vtkXMLUnstructuredGridReader (0x1398ce0): Cannot read point data array "U" from PointData in piece 0. The data array in the element may be tooshort.I noticed that the reader crashes on some absurd numbers appearing inPointData "U" as [...] -3.949617e-105 9.769857e-104 0.000000e+00 [...]If i set those to zero, there is no problem. but the appearance of such 'almost-zero' numbers is clearly connected to using the parameter C.do you suggest me another way how to do this?or do you think that the VTK output writing should check for numbers < DOLFIN_EPS to avoid such corrupted output, and the approach is ok?best regards and thanks for the support, patrick-- Anders_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@xxxxxxxxxxxxxxxxxxx Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@xxxxxxxxxxxxxxxxxxx Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@xxxxxxxxxxxxxxxxxxx Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@xxxxxxxxxxxxxxxxxxx Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@xxxxxxxxxxxxxxxxxxx Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@xxxxxxxxxxxxxxxxxxx Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
| Thread Previous • Date Previous • Date Next • Thread Next |
