Use locks only when modifying entries, not whenever accessing through get_prolongatio...

git-svn-id: https://svn.dealii.org/trunk@29582 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/source/fe/fe_q_base.cc b/deal.II/source/fe/fe_q_base.cc
index 0e106d7c2e9379eeb123db729b2583c3feeacd96..43d46d3dc93da1a4e2ea1858bb82447e393af973 100644 (file)
--- a/deal.II/source/fe/fe_q_base.cc
+++ b/deal.II/source/fe/fe_q_base.cc
@@ -1020,9 +1020,15 @@ FE_Q_Base<POLY,dim,spacedim>
           ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
 
   // initialization upon first request
-  Threads::Mutex::ScopedLock lock(this->mutex);
   if (this->prolongation[refinement_case-1][child].n() == 0)
     {
+      Threads::Mutex::ScopedLock lock(this->mutex);
+
+      // if matrix got updated while waiting for the lock
+      if (this->prolongation[refinement_case-1][child].n() ==
+          this->dofs_per_cell)
+        return this->prolongation[refinement_case-1][child];
+
       // distinguish q/q_dg0 case: only treat Q dofs first
       const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
 
@@ -1072,9 +1078,7 @@ FE_Q_Base<POLY,dim,spacedim>
           }
       }
 
-      FullMatrix<double>& prolongate = const_cast<FullMatrix<double>&>
-        (this->prolongation[refinement_case-1][child]);
-      prolongate.reinit(this->dofs_per_cell, this->dofs_per_cell);
+      FullMatrix<double> prolongate (this->dofs_per_cell, this->dofs_per_cell);
 
       // go through the points in diagonal to capture variation in all
       // directions simultaneously
@@ -1167,6 +1171,10 @@ FE_Q_Base<POLY,dim,spacedim>
           Assert (std::fabs(sum-1.) < eps, ExcInternalError());
         }
 #endif
+
+      // swap matrices
+      std::swap(const_cast<FullMatrix<double> &>
+                (this->prolongation[refinement_case-1][child]), prolongate);
     }
 
   // finally return the matrix
@@ -1189,12 +1197,16 @@ FE_Q_Base<POLY,dim,spacedim>
           ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
 
   // initialization upon first request
-  Threads::Mutex::ScopedLock lock(this->mutex);
   if (this->restriction[refinement_case-1][child].n() == 0)
     {
-      FullMatrix<double> &restriction =
-        const_cast<FullMatrix<double> &>(this->restriction[refinement_case-1][child]);
+      Threads::Mutex::ScopedLock lock(this->mutex);
+
+      // if matrix got updated while waiting for the lock...
+      if (this->restriction[refinement_case-1][child].n() ==
+          this->dofs_per_cell)
+        return this->restriction[refinement_case-1][child];
 
+      FullMatrix<double> restriction(this->dofs_per_cell, this->dofs_per_cell);
       // distinguish q/q_dg0 case
       const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
 
@@ -1291,7 +1303,13 @@ FE_Q_Base<POLY,dim,spacedim>
             restriction(this->dofs_per_cell-1,this->dofs_per_cell-1) =
               1./GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case));
         }
+
+      // swap matrices
+      std::swap(const_cast<FullMatrix<double> &>
+                (this->restriction[refinement_case-1][child]), restriction);
+
     }
+
   return this->restriction[refinement_case-1][child];
 }
 

diff --git a/deal.II/source/fe/fe_system.cc b/deal.II/source/fe/fe_system.cc
index 76c79d2b8ef0e0385af5514e83487953e4e70afb..f20f55657b150c3307f7745a6d9cd5c48be301e4 100644 (file)
--- a/deal.II/source/fe/fe_system.cc
+++ b/deal.II/source/fe/fe_system.cc
@@ -607,6 +607,13 @@ FESystem<dim,spacedim>
   Threads::Mutex::ScopedLock lock(this->mutex);
   if (this->restriction[refinement_case-1][child].n() == 0)
     {
+      Threads::Mutex::ScopedLock lock(this->mutex);
+
+      // check if updated while waiting for lock
+      if (this->restriction[refinement_case-1][child].n() ==
+          this->dofs_per_cell)
+        return this->restriction[refinement_case-1][child];
+
       // Check if some of the matrices of the base elements are void.
       bool do_restriction = true;
 
@@ -627,10 +634,8 @@ FESystem<dim,spacedim>
       // if we did not encounter void matrices, initialize the matrix sizes
       if (do_restriction)
         {
-          FullMatrix<double> &restriction =
-            const_cast<FullMatrix<double>&>(this->restriction[refinement_case-1]
-                                            [child]);
-          restriction.reinit(this->dofs_per_cell, this->dofs_per_cell);
+          FullMatrix<double> restriction(this->dofs_per_cell,
+                                         this->dofs_per_cell);
 
           // distribute the matrices of the base finite elements to the
           // matrices of this object. for this, loop over all degrees of
@@ -664,6 +669,9 @@ FESystem<dim,spacedim>
                 // entries of the matrices:
                 restriction(i,j) = (*base_matrices[base])(base_index_i,base_index_j);
               }
+
+          std::swap(restriction, const_cast<FullMatrix<double> &>
+                    (this->restriction[refinement_case-1][child]));
         }
     }
 
@@ -685,17 +693,19 @@ FESystem<dim,spacedim>
   Assert (child<GeometryInfo<dim>::n_children(refinement_case),
           ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
 
-  // initialization upon first request
-  Threads::Mutex::ScopedLock lock(this->mutex);
+  // initialization upon first request, construction completely analogous to
+  // restriction matrix
   if (this->prolongation[refinement_case-1][child].n() == 0)
     {
-      // Check if some of the matrices of the base elements are void.
-      bool do_prolongation = true;
+      Threads::Mutex::ScopedLock lock(this->mutex);
+
+      if (this->prolongation[refinement_case-1][child].n() ==
+          this->dofs_per_cell)
+        return this->prolongation[refinement_case-1][child];
 
-      // shortcut for accessing local prolongations further down
+      bool do_prolongation = true;
       std::vector<const FullMatrix<double> *>
         base_matrices(this->n_base_elements());
-
       for (unsigned int i=0; i<this->n_base_elements(); ++i)
         {
           base_matrices[i] =
@@ -706,45 +716,27 @@ FESystem<dim,spacedim>
       Assert(do_prolongation,
              (typename FiniteElement<dim,spacedim>::ExcEmbeddingVoid()));
 
-
       if (do_prolongation)
         {
-          FullMatrix<double> &prolongate =
-            const_cast<FullMatrix<double> &>(this->prolongation[refinement_case-1][child]);
-          prolongate.reinit(this->dofs_per_cell, this->dofs_per_cell);
+          FullMatrix<double> prolongate (this->dofs_per_cell,
+                                         this->dofs_per_cell);
 
-          // distribute the matrices of the base finite elements to the
-          // matrices of this object. for this, loop over all degrees of
-          // freedom and take the respective entry of the underlying base
-          // element.
-          //
-          // note that we by definition of a base element, they are
-          // independent, i.e. do not couple. only DoFs that belong to the
-          // same instance of a base element may couple
           for (unsigned int i=0; i<this->dofs_per_cell; ++i)
             for (unsigned int j=0; j<this->dofs_per_cell; ++j)
               {
-                // first find out to which base element indices i and j
-                // belong, and which instance thereof in case the base element
-                // has a multiplicity greater than one. if they should not
-                // happen to belong to the same instance of a base element,
-                // then they cannot couple, so go on with the next index
                 if (this->system_to_base_table[i].first !=
                     this->system_to_base_table[j].first)
                   continue;
-
-                // so get the common base element and the indices therein:
                 const unsigned int
                   base = this->system_to_base_table[i].first.first;
 
                 const unsigned int
                   base_index_i = this->system_to_base_table[i].second,
                   base_index_j = this->system_to_base_table[j].second;
-
-                // if we are sure that DoFs i and j may couple, then copy
-                // entries of the matrices:
                 prolongate(i,j) = (*base_matrices[base])(base_index_i,base_index_j);
               }
+          std::swap(prolongate, const_cast<FullMatrix<double> &>
+                    (this->prolongation[refinement_case-1][child]));
         }
     }