]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement the building of interface constraints for system-fe's also.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 12 Apr 1999 20:06:28 +0000 (20:06 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 12 Apr 1999 20:06:28 +0000 (20:06 +0000)
git-svn-id: https://svn.dealii.org/trunk@1129 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_system.h
deal.II/deal.II/source/fe/fe_system.cc

index 9afcc3b28f1be57c2d59278d3304457a0722870d..073ef326b7c883c1b84745ece9467f464d036eda 100644 (file)
@@ -474,6 +474,11 @@ class FESystem //<dim>
                                      * Used by #initialize#.
                                      */
     void build_face_table();
+
+                                    /**
+                                     * Used by #initialize#.
+                                     */
+    void build_interface_constraints ();
     
                                     /**
                                      *Exception.
index 329b252dc8c782560a4e671f261ded866b36cc18..04bf948dc25953834d586c9c0a072dbba7856bcb 100644 (file)
@@ -256,6 +256,120 @@ FESystem<dim>::build_face_table()
 
 
 
+template <int dim>
+void FESystem<dim>::build_interface_constraints () 
+{
+                                  // the layout of the constraints matrix is
+                                  // described in the FiniteElement class. you
+                                  // may want to look there first before trying
+                                  // to understand the following, especially
+                                  // the mapping of the #n# index.
+                                  //
+                                  // in order to map it to the fe-system class,
+                                  // we have to know which base element a
+                                  // degree of freedom within a vertex, line,
+                                  // etc belongs to. this can be accomplished
+                                  // by the system_to_component_index
+                                  // function in conjunction with the
+                                  // numbers first_{line,quad,...}_index
+  for (unsigned int n=0; n<interface_constraints.n(); ++n)
+    for (unsigned int m=0; m<interface_constraints.m(); ++m)
+      {
+                                        // for the pair (n,m) find out which
+                                        // component they belong to and
+                                        // the number therein
+                                        //
+                                        // first value in pair is component,
+                                        // second is index
+       const pair<unsigned int, unsigned int> n_index
+         = face_system_to_component_index (n);
+
+       pair<unsigned int, unsigned int> m_index;
+       switch (dim)
+         {
+           case 1:
+           {
+                                              // we should never get here!
+                                              // (in 1d, the constraints matrix
+                                              // should be of size zero)
+             Assert (false, ExcInternalError());
+             break;
+           };
+
+           case 2:
+           {
+                                              // the indices m=0..d_v-1 are
+                                              // from the center vertex.
+                                              // their order is the same
+                                              // as for the first vertex
+                                              // of the whole cell, so we
+                                              // can use the
+                                              // system_to_component_index
+                                              // function (using the
+                                              // face_s_t_c_i function would
+                                              // yield the same)
+             if (m < dofs_per_vertex)
+               m_index = system_to_component_index (m);
+             else
+                                                // then come the two sets of
+                                                // line indices
+               {
+                 const unsigned int index_in_line
+                   = (m-dofs_per_vertex) % dofs_per_line;
+                 const unsigned int sub_line
+                   = (m-dofs_per_vertex) / dofs_per_line;
+                 Assert (sub_line < 2, ExcInternalError());
+                 
+                                                  // get the component by
+                                                  // asking s_t_c_index and
+                                                  // tweaking the index a bit
+                 m_index.first = system_to_component_index
+                                 (GeometryInfo<2>::vertices_per_cell * dofs_per_vertex
+                                  + index_in_line).first;
+                                                  // first find out the how-many'th
+                                                  // line index of that component
+                                                  // this was
+                 m_index.second = (system_to_component_index
+                                   (GeometryInfo<2>::vertices_per_cell * dofs_per_vertex
+                                    + index_in_line).second
+                                   - base_element (m_index.first).first_line_index)
+                                                  // then add the number of dofs
+                                                  // per vertex to get the index
+                                                  // on the first line
+                                  + base_element(m_index.first).dofs_per_vertex
+                                                  // if on the second line: add
+                                                  // some more
+                                  + base_element(m_index.first).dofs_per_line * sub_line;
+               };
+             break;
+           };
+
+           case 3:
+                                                  // it should be straightforward
+                                                  // to implement this for 3d from
+                                                  // the code above, but I haven't
+                                                  // done it yet.
+                 
+           default:
+                 Assert (false, ExcNotImplemented());
+         };
+
+                                        // now that we gathered all information:
+                                        // use it to build the matrix. note
+                                        // that if n and m belong to different
+                                        // components, there definitely will be
+                                        // no coupling
+       if (n_index.first == m_index.first)
+         interface_constraints(m,n)
+           = base_element(n_index.first).constraints()(m_index.second,
+                                                       n_index.second);
+      };
+};
+
+
+
+
+
 template <int dim>
 void FESystem<dim>::initialize ()
 {
@@ -265,26 +379,27 @@ void FESystem<dim>::initialize ()
                                   // finite elements to the matrices of
                                   // this object
   for (unsigned int component=0; component<n_components; ++component)
+                                    // transform restriction and
+                                    // prolongation matrices
     for (unsigned int i=0; i<base_element(component_to_base_table[component]).total_dofs; ++i)
       for (unsigned int j=0; j<base_element(component_to_base_table[component]).total_dofs; ++j)
                                         // only fill block diagonals, no
                                         // intermixing of subelements
-       {
-         for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
-           {
-             restriction[child] (component_to_system_index (component,i),
-                                 component_to_system_index (component, j))
-               = base_element(component_to_base_table[component]).restrict(child)(i,j);
-             prolongation[child] (component_to_system_index (component,i),
-                                  component_to_system_index (component, j))
-               = base_element(component_to_base_table[component]).prolongate(child)(i,j);
-           };
-
-// still to do:          
-//       interface_constraints (// component_to_system_index (component,i),
-//                              component_to_system_index (component, j))
-//         = base_element(component_to_base_table[component]).constraints()(i,j);
-       };
+       for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
+         {
+           restriction[child] (component_to_system_index (component,i),
+                               component_to_system_index (component, j))
+             = base_element(component_to_base_table[component]).restrict(child)(i,j);
+           prolongation[child] (component_to_system_index (component,i),
+                                component_to_system_index (component, j))
+             = base_element(component_to_base_table[component]).prolongate(child)(i,j);
+         };
+
+
+                                  // now set up the interface constraints.
+                                  // this is kind'o hairy, so don't try
+                                  // to do it dimension independent
+  build_interface_constraints ();
 };
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.