#10 open
Jed Brown

Generic boundary conditions

Reported by Jed Brown | December 1st, 2008 @ 03:43 PM

We need a generic way to impose a wide range of boundary conditions including inhomogeneous Dirichlet (e.g. no-slip), inhomogeneous Neumann (e.g. stress), slip (nonlinear Robin in tangent plane, inhomogeneous Dirichlet in normal direction).

My proposal

  • User creates a scalar valued function space for the entire domain

  • User obtains array of entities for each boundary type (usually by looking up tags on the mesh), with per-face tag indicating whether normal is "in" or "out".

  • User registers set, name, and boundary constraint function with function space. The boundary constraint function has a prototype something like


dErr UserBdyConstrain(void *ctx,const dReal x[3],const dReal basis[3][3],dInt *nglobal,dReal T[]);

If the FS has D dofs per node, T should be a D by D matrix which writes the new dofs as a linear combination of the old dofs. In all cases I can think of, the matrix T should be orthogonal (the implementation should probably confirm this). The first nglobal rows will be global dofs and the rest will be Dirichlet values. The Dirichlet values may be set in function evaluation. As an example, suppose we are building a combined velocity/temperature space with a slip boundary.


const dReal *n = basis[0],*t = basis[1],*tt = basis[2];
const dReal mat[4][4] = {{t[0],t[1],t[2],0},{tt[0],tt[1],tt[2],0},{0,0,0,1},{n[0],n[1],n[2],0}};
memcpy(T,mat,sizeof(mat));
*nglobal = 3;

The implementation will represent the first three values in the global basis (tangent velocity and temperature) and one component (normal velocity) Dirichlet. The implementation could choose not to store entries in the matrix which are exactly zero. The user can get a vector of all Dirichlet values from the FS and scatter to/from it in function evaluation. Except by scattering into the standard basis, this Dirichlet vector will be opaque to the user.

  • During function evaluation and Jacobian application, the user needs to iterate over the faces in each boundary type and do element operations. For this, the FS will extend the expanded vector to include chunks for each boundary condition and create dRule and dEFS objects for each non-Dirichlet face. Then, much like dFSGetElements, the user will call dFSGetFaces and compute the boundary integrals. The FS will provide face geometry and normal vectors.

Comments and changes to this ticket

  • Jed Brown

    Jed Brown December 1st, 2008 @ 03:50 PM

    • State changed from “new” to “open”
  • Jed Brown

    Jed Brown March 26th, 2009 @ 08:51 PM

    Upon further thought, the design proposed above is somewhat limiting and confusing. While I still agree with the constraint matrix part, we need a way of handling association with geometric entities.

    ITAPS boundary conditions

    Boundary conditions need to associate mesh entities with geometric entities (which provide a high-accuracy representation with normal and tangent vectors, and curvature if necessary). Cubit boundaries are represented in MOAB by tagging sets with an integer id NEUMANN_SET. A tagged set with boundary condition 100, say, will contain multiple subsets. There might be one subset with the tag SENSE and value -1. With the exception of the sense set, each subset of NEUMANN_SET=100 contains mesh entities and an association with a geometric entity. If it is also in NEUMANN_SET=100:SENSE=-1, then the normals of the faces point inwards instead of outwards.

    For difficult boundary conditions, such as for implicit slip conditions, it is useful for user code (function evaluation) to have access to the geometric associations. This suggests that there is benefit to tagging the mesh entities with necessary metadata. Note that some of this metadata cannot be exchanged in parallel, it would only make sense on the PE it is created on.

    Matrix performance and block formats (BAIJ)

    There is a performance benefit to using the BAIJ matrix format. While the use of inodes nullifies the computation/memory improvement for matrix multiplication and factorization, the ability to fully unroll these operations does make a performance difference. In addition MatPBRelax is a stronger smoother than MatRelax and there is no good analogue for non-block formats.

    Imposing Dirichlet (essential) conditions with BAIJ

    To impose Dirichlet conditions on all components in a block, the entire block should be removed from the global system, the Dirichlet values are stored separately, and function evaluation will add these together in the expanded space. To impose Dirichlet conditions on only some components of the block, the constrained rows and columns of the Jacobian should be zero. During assembly, this can be achieved by zeroing them on matrix insertion. For function evaluation, this means scattering in correct values before local operations, and setting those components of the residual vector to the difference between the current iterate and the correct values.

    Imposing weak conditions

    Homogeneous Neumann conditions are the "do nothing" condition. Inhomogeneous conditions require face integrals to be performed. If the condition is of Robin type, is a "no condition" outflow, or is nonlinear, it will contribute to the Jacobian. In the latter case, there will be an additional term arising from differentiation.

    Proposal

    The user registers boundaries, indicating constraint information (all components constrained, none constrained, some constrained) and whether the faces should be represented in the expanded space (i.e. surface integrals must be evaluated).

    Storage

    An extension to VecGhost: one large array holds [owned, ghosted, dirichlet] blocks. The global form is just the first bit, the local form includes ghosted values (and is a sequential vector) and the full form is identical to the local form but includes Dirichlet values. A separate Dirichlet vector is allocated and a scatter to the full form is defined.

    Setting up the function space

    Dohp tags entities with the entity local index relative to this function space. User code will have access to the iMesh sets, hence it will be able to obtain geometry information if necessary. This index can be used to find the starting offset in the expanded vector (for evaluating weak forms) and the starting offset in the full vector (for strong manipulation). The index is also used to get the dEFS needed for evaluating weak forms.

  • Jed Brown

    Jed Brown April 17th, 2009 @ 04:27 PM

    As an addendum, the proposed storage scheme doesn't work well since it is desirable to have a global form of the Dirichlet vector.

    A different proposal (mostly implemented)

    The global closure vector is [owned, owned Dirichlet]. There are both global and global Dirichlet vectors (the latter is useful to ensure consistency in parallel). These are both ghosted vectors and there are scatters from the global closure to each global form.

  • Jed Brown

    Jed Brown April 17th, 2009 @ 08:30 PM

    (from [69efc49fcd9c59a7eff1850a61b24f4b7eb7247a]) Work on handling boundary conditions. [#5] [#10]

    • Provide Global, Dirichlet, and Closure vectors. Global and Dirichlet vectors are VecGhost so they have a local form. Scatters are provided from Closure to the others.

    • Use MAIJ format for element assembly matrices (allows vector-valued spaces)

    • Keep track of boundaries by using ITAPS/MOAB tag conventions instead of private data structure (this is useful since the user needs to see the set structure if they need to relate to geometric entities).

    • Add dFSHomogeneousMode to handle which boundary values should be mapped into expanded space. http://github.com/jedbrown/dohp/...

Please Sign in or create a free account to add a new ticket.

With your very own profile, you can contribute to projects, track your activity, watch tickets, receive and update tickets through your email and much more.

New-ticket Create new ticket

Create your profile

Help contribute to this project by taking a few moments to create your personal profile. Create your profile ยป

An implementation of the ``dual order hp'' version of the finite element method. This project targets parallel domain-decomposition methods for strongly coupled nonlinear problems with PDE constraints.

People watching this ticket

Tags

Referenced by

Pages