← Back to branch summary

~ubuntu-branches/ubuntu/vivid/deal.ii/vivid

« back to all changes in this revision

Viewing changes to deal.II/source/fe/fe_tools.cc

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV, Adam C. Powell, IV, Denis Barbier
  • Date: 2010-07-29 13:47:01 UTC
  • mfrom: (3.1.3 sid)
  • Revision ID: james.westby@ubuntu.com-20100729134701-akb8jb3stwge8tcm
Tags: 6.3.1-1
http://bugs.debian.org/562332
[ Adam C. Powell, IV ]
* Changed to source format 3.0 (quilt).
* Changed maintainer to debian-science with Adam Powell as uploader.
* Added source lintian overrides about Adam Powell's name.
* Added Vcs info on git repository.
* Bumped Standards-Version.
* Changed stamp-patch to patch target and fixed its application criterion.
* Moved make_dependencies and expand_instantiations to a versioned directory
  to avoid shlib package conflicts.

[ Denis Barbier ]
* New upstream release (closes: #562332).
  + Added libtbb support.
  + Forward-ported all patches.
* Updates for new PETSc version, including workaround for different versions
  of petsc and slepc.
* Add debian/watch.
* Update to debhelper 7.
* Added pdebuild patch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
deal.II/source/fe/fe_tools.cc
1
1
//---------------------------------------------------------------------------
2
 
//    $Id: fe_tools.cc 18724 2009-04-23 23:06:37Z bangerth $
 
2
//    $Id: fe_tools.cc 20944 2010-04-06 04:13:47Z bangerth $
3
3
//    Version: $Name$
4
4
//
5
 
//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
 
5
//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
6
6
//
7
7
//    This file is subject to QPL and may not be  distributed
8
8
//    without copyright and license information. Please refer
51
51
DEAL_II_NAMESPACE_OPEN
52
52
 
53
53
 
54
 
namespace 
 
54
namespace
55
55
{
56
56
                                   // a shared pointer to factory
57
57
                                   // objects, to ensure that they get
62
62
  typedef
63
63
  std_cxx1x::shared_ptr<const FETools::FEFactoryBase<deal_II_dimension> >
64
64
  FEFactoryPointer;
65
 
  
 
65
 
66
66
                                   // a function that returns the
67
67
                                   // default set of finite element
68
68
                                   // names and factory objects for
70
70
                                   // fe_name_map below
71
71
#ifdef DEAL_II_ANON_NAMESPACE_BUG
72
72
    static
73
 
#endif  
 
73
#endif
74
74
  std::map<std::string,FEFactoryPointer>
75
75
  get_default_fe_names ()
76
76
  {
101
101
 
102
102
    return default_map;
103
103
  }
104
 
  
105
 
 
106
 
  
 
104
 
 
105
 
 
106
 
107
107
                                   // have a lock that guarantees that
108
108
                                   // at most one thread is changing
109
109
                                   // and accessing the fe_name_map
152
152
 
153
153
 
154
154
 
155
 
namespace 
 
155
namespace
156
156
{
157
157
 
158
158
                                   // forwarder function for
176
176
    fe2.get_interpolation_matrix (fe1, interpolation_matrix);
177
177
  }
178
178
 
179
 
  
 
179
 
180
180
  template <int dim, typename number, int spacedim>
181
181
  inline
182
182
  void gim_forwarder (const FiniteElement<dim,spacedim> &fe1,
221
221
 
222
222
    Assert (dim<10, ExcNotImplemented());
223
223
    const char dim_char = '0'+dim;
224
 
    
 
224
 
225
225
    if ((position+3 < name.size())
226
226
        &&
227
227
        (name[position] == '<')
252
252
  Assert(renumbering.size() == element.dofs_per_cell,
253
253
         ExcDimensionMismatch(renumbering.size(),
254
254
                              element.dofs_per_cell));
255
 
  
 
255
 
256
256
  comp_start.resize(element.n_base_elements());
257
 
  
 
257
 
258
258
  unsigned int k=0;
259
259
  for (unsigned int i=0;i<comp_start.size();++i)
260
260
    {
261
261
      comp_start[i].resize(element.element_multiplicity(i));
262
262
      const unsigned int increment
263
263
        = element.base_element(i).dofs_per_cell;
264
 
      
 
264
 
265
265
      for (unsigned int j=0;j<comp_start[i].size();++j)
266
266
        {
267
267
          comp_start[i][j] = k;
268
268
          k += increment;
269
269
        }
270
270
    }
271
 
  
 
271
 
272
272
                                   // For each index i of the
273
273
                                   // unstructured cellwise
274
274
                                   // numbering, renumbering
289
289
void FETools::compute_block_renumbering (
290
290
  const FiniteElement<dim,spacedim>& element,
291
291
  std::vector<unsigned int>& renumbering,
292
 
  std::vector<unsigned int>& start_indices)
 
292
  std::vector<unsigned int>& block_data,
 
293
  bool return_start_indices)
293
294
{
294
295
  Assert(renumbering.size() == element.dofs_per_cell,
295
296
         ExcDimensionMismatch(renumbering.size(),
296
297
                              element.dofs_per_cell));
297
 
  Assert(start_indices.size() == element.n_blocks(),
298
 
         ExcDimensionMismatch(start_indices.size(),
 
298
  Assert(block_data.size() == element.n_blocks(),
 
299
         ExcDimensionMismatch(block_data.size(),
299
300
                              element.n_blocks()));
300
 
  
 
301
 
301
302
  unsigned int k=0;
302
303
  unsigned int i=0;
303
304
  for (unsigned int b=0;b<element.n_base_elements();++b)
304
305
    for (unsigned int m=0;m<element.element_multiplicity(b);++m)
305
306
      {
306
 
        start_indices[i++] = k;
 
307
        block_data[i++] = (return_start_indices)
 
308
                             ? k
 
309
                             : (element.base_element(b).n_dofs_per_cell());
307
310
        k += element.base_element(b).n_dofs_per_cell();
308
311
      }
309
312
  Assert (i == element.n_blocks(), ExcInternalError());
310
 
//TODO:[GK] This does not work for a single RT  
 
313
 
 
314
  std::vector<unsigned int> start_indices(block_data.size());
 
315
  k = 0;
 
316
  for (unsigned int i=0;i<block_data.size();++i)
 
317
    if (return_start_indices)
 
318
      start_indices[i] = block_data[i];
 
319
    else
 
320
      {
 
321
        start_indices[i] = k;
 
322
        k += block_data[i];
 
323
      }
 
324
 
 
325
//TODO:[GK] This does not work for a single RT
311
326
  for (unsigned int i=0;i<element.dofs_per_cell;++i)
312
327
    {
313
328
      std::pair<unsigned int, unsigned int>
337
352
                                   // the FE wants to implement things
338
353
                                   // itself:
339
354
  bool fe_implements_interpolation = true;
340
 
  try 
 
355
  try
341
356
    {
342
357
      gim_forwarder (fe1, fe2, interpolation_matrix);
343
358
    }
367
382
  Assert(fe2_support_points.size()==fe2.dofs_per_cell,
368
383
         typename FEL::ExcFEHasNoSupportPoints());
369
384
 
370
 
  for (unsigned int i=0; i<fe2.dofs_per_cell; ++i)      
 
385
  for (unsigned int i=0; i<fe2.dofs_per_cell; ++i)
371
386
    {
372
387
      const unsigned int i1 = fe2.system_to_component_index(i).first;
373
388
      for (unsigned int j=0; j<fe1.dofs_per_cell; ++j)
377
392
            interpolation_matrix(i,j) = fe1.shape_value (j,fe2_support_points[i]);
378
393
          else
379
394
            interpolation_matrix(i,j)=0.;
380
 
        }  
 
395
        }
381
396
    }
382
397
}
383
398
 
391
406
  Assert (fe1.n_components() == fe2.n_components(),
392
407
          ExcDimensionMismatch(fe1.n_components(), fe2.n_components()));
393
408
  Assert(interpolation_matrix.m()==fe1.dofs_per_cell &&
394
 
         interpolation_matrix.n()==fe1.dofs_per_cell, 
 
409
         interpolation_matrix.n()==fe1.dofs_per_cell,
395
410
         ExcMatrixDimensionMismatch(interpolation_matrix.m(),
396
411
                                    interpolation_matrix.n(),
397
412
                                    fe1.dofs_per_cell,
398
413
                                    fe1.dofs_per_cell));
399
 
    
 
414
 
400
415
  FullMatrix<number> first_matrix (fe2.dofs_per_cell, fe1.dofs_per_cell);
401
416
  FullMatrix<number> second_matrix(fe1.dofs_per_cell, fe2.dofs_per_cell);
402
 
  
 
417
 
403
418
  get_interpolation_matrix(fe1, fe2, first_matrix);
404
419
  get_interpolation_matrix(fe2, fe1, second_matrix);
405
420
 
417
432
  Assert (fe1.n_components() == fe2.n_components(),
418
433
          ExcDimensionMismatch(fe1.n_components(), fe2.n_components()));
419
434
  Assert(difference_matrix.m()==fe1.dofs_per_cell &&
420
 
         difference_matrix.n()==fe1.dofs_per_cell, 
 
435
         difference_matrix.n()==fe1.dofs_per_cell,
421
436
         ExcMatrixDimensionMismatch(difference_matrix.m(),
422
437
                                    difference_matrix.n(),
423
438
                                    fe1.dofs_per_cell,
424
439
                                    fe1.dofs_per_cell));
425
 
   
 
440
 
426
441
  FullMatrix<number> interpolation_matrix(fe1.dofs_per_cell);
427
442
  get_back_interpolation_matrix(fe1, fe2, interpolation_matrix);
428
 
  
 
443
 
429
444
  for (unsigned int i=0; i<fe1.dofs_per_cell; ++i)
430
445
    difference_matrix(i,i) = 1.;
431
 
  
 
446
 
432
447
                                   // compute difference
433
448
  difference_matrix.add (-1, interpolation_matrix);
434
449
}
448
463
                                    fe2.dofs_per_cell,
449
464
                                    fe1.dofs_per_cell));
450
465
  matrix = 0;
451
 
  
 
466
 
452
467
  unsigned int n1 = fe1.dofs_per_cell;
453
468
  unsigned int n2 = fe2.dofs_per_cell;
454
469
 
456
471
                                   // the unit cell
457
472
  Triangulation<dim,spacedim> tr;
458
473
  GridGenerator::hyper_cube(tr);
459
 
  
 
474
 
460
475
                                   // Choose a quadrature rule
461
476
                                   // Gauss is exact up to degree 2n-1
462
477
  const unsigned int degree = std::max(fe1.tensor_degree(), fe2.tensor_degree());
463
478
  Assert (degree != numbers::invalid_unsigned_int,
464
479
          ExcNotImplemented());
465
 
  
 
480
 
466
481
  QGauss<dim> quadrature(degree+1);
467
482
                                   // Set up FEValues.
468
483
  const UpdateFlags flags = update_values | update_quadrature_points | update_JxW_values;
490
505
                                   // mass matrix to be
491
506
                                   // well-conditioned
492
507
  mass.gauss_jordan();
493
 
  
 
508
 
494
509
                                   // Now, test every function of fe1
495
510
                                   // with test functions of fe2 and
496
511
                                   // compute the projection of each
497
512
                                   // unit vector.
498
513
  Vector<double> b(n2);
499
514
  Vector<double> x(n2);
500
 
  
 
515
 
501
516
  for (unsigned int j=0;j<n1;++j)
502
517
    {
503
518
      b = 0.;
509
524
            const double v = val2.shape_value(i,k);
510
525
            b(i) += u*v*w;
511
526
          }
512
 
      
 
527
 
513
528
                                       // Multiply by the inverse
514
529
      mass.vmult(x,b);
515
530
      for (unsigned int i=0;i<n2;++i)
530
545
  Assert (N.m()==n_dofs, ExcDimensionMismatch(N.m(), n_dofs));
531
546
 
532
547
  const std::vector<Point<dim> >& points = fe.get_generalized_support_points();
533
 
  
 
548
 
534
549
                                   // We need the values of the
535
550
                                   // polynomials in all generalized
536
551
                                   // support points.
537
552
  std::vector<std::vector<double> >
538
553
    values (dim, std::vector<double>(points.size()));
539
 
  
 
554
 
540
555
                                   // In this vector, we store the
541
556
                                   // result of the interpolation
542
557
  std::vector<double> local_dofs(n_dofs);
543
 
  
 
558
 
544
559
                                   // One row per shape
545
560
                                   // function. Remember that these
546
561
                                   // are the 'raw' shape functions
598
613
  unsigned int ref_case = (isotropic_only)
599
614
                          ? RefinementCase<dim>::isotropic_refinement
600
615
                          : RefinementCase<dim>::cut_x;
601
 
  
 
616
 
602
617
  for (;ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
603
618
    {
604
619
      const unsigned int nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
607
622
          Assert(matrices[ref_case-1][i].n() == n, ExcDimensionMismatch(matrices[ref_case-1][i].n(),n));
608
623
          Assert(matrices[ref_case-1][i].m() == n, ExcDimensionMismatch(matrices[ref_case-1][i].m(),n));
609
624
        }
610
 
  
 
625
 
611
626
                                   // Set up meshes, one with a single
612
627
                                   // reference cell and refine it once
613
628
      Triangulation<dim,spacedim> tria;
618
633
      MappingCartesian<dim> mapping;
619
634
      QGauss<dim> q_fine(degree+1);
620
635
      const unsigned int nq = q_fine.size();
621
 
  
 
636
 
622
637
      FEValues<dim> fine (mapping, fe, q_fine,
623
638
                          update_quadrature_points | update_JxW_values | update_values);
624
 
  
 
639
 
625
640
                                       // We search for the polynomial on
626
641
                                       // the small cell, being equal to
627
642
                                       // the coarse polynomial in all
628
643
                                       // quadrature points.
629
 
                                    
 
644
 
630
645
                                       // First build the matrix for this
631
646
                                       // least squares problem. This
632
647
                                       // contains the values of the fine
633
648
                                       // cell polynomials in the fine
634
649
                                       // cell grid points.
635
 
                                    
 
650
 
636
651
                                       // This matrix is the same for all
637
652
                                       // children.
638
653
      fine.reinit(tria.begin_active());
639
654
      FullMatrix<number> A(nq*nd, n);
640
 
      for (unsigned int d=0;d<nd;++d)
641
 
        for (unsigned int k=0;k<nq;++k)
642
 
          for (unsigned int j=0;j<n;++j)
 
655
      for (unsigned int j=0;j<n;++j)
 
656
        for (unsigned int d=0;d<nd;++d)
 
657
          for (unsigned int k=0;k<nq;++k)
643
658
            A(k*nd+d,j) = fine.shape_value_component(j,k,d);
644
659
 
645
660
      Householder<double> H(A);
646
 
  
 
661
 
647
662
      Vector<number> v_coarse(nq*nd);
648
663
      Vector<number> v_fine(n);
649
 
  
 
664
 
650
665
      unsigned int cell_number = 0;
651
666
      for (typename Triangulation<dim>::active_cell_iterator fine_cell
652
667
             = tria.begin_active();
664
679
          coarse.reinit(tria.begin(0));
665
680
 
666
681
          FullMatrix<double> &this_matrix = matrices[ref_case-1][cell_number];
667
 
      
 
682
          v_coarse = 0;
 
683
 
668
684
                                           // Compute this once for each
669
685
                                           // coarse grid basis function
670
686
          for (unsigned int i=0;i<n;++i)
675
691
                                               // function values of the
676
692
                                               // coarse grid function in
677
693
                                               // each quadrature point.
678
 
              for (unsigned int d=0;d<nd;++d)
679
 
                for (unsigned int k=0;k<nq;++k)
680
 
                  v_coarse(k*nd+d) = coarse.shape_value_component(i,k,d);
 
694
              if (fe.is_primitive())
 
695
                {
 
696
                  const unsigned int d = fe.system_to_component_index(i).first;
 
697
                  const double * phi_i = &coarse.shape_value (i,0);
 
698
                  for (unsigned int k=0;k<nq;++k)
 
699
                    v_coarse(k*nd+d) = phi_i[k];
 
700
                }
 
701
              else
 
702
                for (unsigned int d=0;d<nd;++d)
 
703
                  for (unsigned int k=0;k<nq;++k)
 
704
                    v_coarse(k*nd+d) = coarse.shape_value_component(i,k,d);
681
705
 
682
706
                                               // solve the least squares
683
707
                                               // problem.
684
708
              const double result = H.least_squares(v_fine, v_coarse);
685
709
              Assert (result < 1.e-12, ExcLeastSquaresError(result));
686
 
            
 
710
 
687
711
                                               // Copy into the result
688
712
                                               // matrix. Since the matrix
689
713
                                               // maps a coarse grid
720
744
  const unsigned int n  = fe.dofs_per_face;
721
745
  const unsigned int nd = fe.n_components();
722
746
  const unsigned int degree = fe.degree;
723
 
  
 
747
 
724
748
  for (unsigned int i=0;i<nc;++i)
725
749
    {
726
750
      Assert(matrices[i].n() == n, ExcDimensionMismatch(matrices[i].n(),n));
727
751
      Assert(matrices[i].m() == n, ExcDimensionMismatch(matrices[i].m(),n));
728
752
    }
729
 
  
 
753
 
730
754
                                   // Set up meshes, one with a single
731
755
                                   // reference cell and refine it once
732
756
  Triangulation<dim,spacedim> tria;
736
760
  MappingCartesian<dim> mapping;
737
761
  QGauss<dim-1> q_gauss(degree+1);
738
762
  const Quadrature<dim> q_fine = QProjector<dim>::project_to_face(q_gauss, face_fine);
739
 
  
 
763
 
740
764
  const unsigned int nq = q_fine.size();
741
765
 
742
766
                                   // In order to make the loops below
791
815
        }
792
816
    }
793
817
  Assert (k == fe.dofs_per_face, ExcInternalError());
794
 
  
 
818
 
795
819
  FEValues<dim> fine (mapping, fe, q_fine,
796
820
                      update_quadrature_points | update_JxW_values | update_values);
797
 
  
 
821
 
798
822
                                   // We search for the polynomial on
799
823
                                   // the small cell, being equal to
800
824
                                   // the coarse polynomial in all
801
825
                                   // quadrature points.
802
 
                                    
 
826
 
803
827
                                   // First build the matrix for this
804
828
                                   // least squares problem. This
805
829
                                   // contains the values of the fine
806
830
                                   // cell polynomials in the fine
807
831
                                   // cell grid points.
808
 
                                    
 
832
 
809
833
                                   // This matrix is the same for all
810
834
                                   // children.
811
835
  fine.reinit(tria.begin_active());
812
836
  FullMatrix<number> A(nq*nd, n);
813
 
  for (unsigned int d=0;d<nd;++d)
814
 
    for (unsigned int k=0;k<nq;++k)
815
 
      for (unsigned int j=0;j<n;++j)
 
837
  for (unsigned int j=0;j<n;++j)
 
838
    for (unsigned int d=0;d<nd;++d)
 
839
      for (unsigned int k=0;k<nq;++k)
816
840
        A(k*nd+d,j) = fine.shape_value_component(face_f_dofs[j],k,d);
817
841
 
818
842
  Householder<double> H(A);
819
 
  
 
843
 
820
844
  Vector<number> v_coarse(nq*nd);
821
845
  Vector<number> v_fine(n);
822
 
  
823
 
  
824
 
  
 
846
 
 
847
 
 
848
 
825
849
  for (unsigned int cell_number = 0; cell_number < GeometryInfo<dim>::max_children_per_face;
826
850
       ++cell_number)
827
851
    {
828
852
      const Quadrature<dim> q_coarse
829
853
        = QProjector<dim>::project_to_subface(q_gauss, face_coarse, cell_number);
830
854
      FEValues<dim> coarse (mapping, fe, q_coarse, update_values);
831
 
      
 
855
 
832
856
      typename Triangulation<dim,spacedim>::active_cell_iterator fine_cell
833
857
        = tria.begin(0)->child(GeometryInfo<dim>::child_cell_on_face(
834
858
          tria.begin(0)->refinement_case(), face_coarse, cell_number));
835
859
      fine.reinit(fine_cell);
836
860
      coarse.reinit(tria.begin(0));
837
 
      
 
861
 
838
862
      FullMatrix<double> &this_matrix = matrices[cell_number];
839
 
      
 
863
 
840
864
                                       // Compute this once for each
841
865
                                       // coarse grid basis function
842
866
      for (unsigned int i=0;i<n;++i)
849
873
                                           // each quadrature point.
850
874
          for (unsigned int d=0;d<nd;++d)
851
875
            for (unsigned int k=0;k<nq;++k)
852
 
              v_coarse(k*nd+d) = coarse.shape_value_component(face_c_dofs[i],k,d);
 
876
              v_coarse(k*nd+d) = coarse.shape_value_component (face_c_dofs[i],k,d);
853
877
 
854
878
                                           // solve the least squares
855
879
                                           // problem.
856
880
          const double result = H.least_squares(v_fine, v_coarse);
857
881
          Assert (result < 1.e-12, ExcLeastSquaresError(result));
858
 
            
 
882
 
859
883
                                           // Copy into the result
860
884
                                           // matrix. Since the matrix
861
885
                                           // maps a coarse grid
904
928
  const unsigned int nd = fe.n_components();
905
929
  const unsigned int degree = fe.degree;
906
930
 
 
931
                                   // prepare FEValues, quadrature etc on
 
932
                                   // coarse cell
 
933
  MappingCartesian<dim> mapping;
 
934
  QGauss<dim> q_fine(degree+1);
 
935
  const unsigned int nq = q_fine.size();
 
936
 
 
937
                                   // create mass matrix on coarse cell.
 
938
  FullMatrix<number> mass(n, n);
 
939
  {
 
940
                                   // set up a triangulation for coarse cell
 
941
    Triangulation<dim,spacedim> tr;
 
942
    GridGenerator::hyper_cube (tr, 0, 1);
 
943
 
 
944
    FEValues<dim> coarse (mapping, fe, q_fine,
 
945
                          update_JxW_values | update_values);
 
946
 
 
947
    typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
 
948
      = tr.begin(0);
 
949
    coarse.reinit (coarse_cell);
 
950
 
 
951
    const std::vector<double> & JxW = coarse.get_JxW_values();
 
952
    for (unsigned int i=0;i<n;++i)
 
953
      for (unsigned int j=0;j<n;++j)
 
954
        if (fe.is_primitive())
 
955
          {
 
956
            const double * coarse_i = &coarse.shape_value(i,0);
 
957
            const double * coarse_j = &coarse.shape_value(j,0);
 
958
            double mass_ij = 0;
 
959
            for (unsigned int k=0;k<nq;++k)
 
960
              mass_ij += JxW[k] * coarse_i[k] * coarse_j[k];
 
961
            mass(i,j) = mass_ij;
 
962
          }
 
963
        else
 
964
          {
 
965
            double mass_ij = 0;
 
966
            for (unsigned int d=0;d<nd;++d)
 
967
              for (unsigned int k=0;k<nq;++k)
 
968
                mass_ij += JxW[k] * coarse.shape_value_component(i,k,d)
 
969
                                  * coarse.shape_value_component(j,k,d);
 
970
            mass(i,j) = mass_ij;
 
971
          }
 
972
 
 
973
                                   // invert mass matrix
 
974
    mass.gauss_jordan();
 
975
  }
 
976
 
907
977
                                   // loop over all possible refinement cases
908
978
  for (unsigned int ref_case = RefinementCase<dim>::cut_x;
909
 
       ref_case<RefinementCase<dim>::isotropic_refinement+1; ++ref_case)
 
979
       ref_case < RefinementCase<dim>::isotropic_refinement+1; ++ref_case)
910
980
    {
911
981
      const unsigned int
912
982
        nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
918
988
          Assert(matrices[ref_case-1][i].m() == n,
919
989
                 ExcDimensionMismatch(matrices[ref_case-1][i].m(),n));
920
990
        }
921
 
  
 
991
 
 
992
                                   // create a respective refinement on the
 
993
                                   // triangulation
922
994
      Triangulation<dim,spacedim> tr;
923
995
      GridGenerator::hyper_cube (tr, 0, 1);
924
996
      tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
925
997
      tr.execute_coarsening_and_refinement();
926
 
  
927
 
      MappingCartesian<dim> mapping;
928
 
      QGauss<dim> q_fine(degree+1);
929
 
      const unsigned int nq = q_fine.size();
930
 
  
931
 
      FEValues<dim> coarse (mapping, fe, q_fine,
932
 
                            update_quadrature_points |
933
 
                            update_JxW_values |
934
 
                            update_values);
 
998
 
935
999
      FEValues<dim> fine (mapping, fe, q_fine,
936
 
                          update_quadrature_points |
937
 
                          update_JxW_values |
 
1000
                          update_quadrature_points | update_JxW_values |
938
1001
                          update_values);
939
 
  
 
1002
 
940
1003
      typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
941
1004
        = tr.begin(0);
942
 
                                       // Compute the coarse level mass
943
 
                                       // matrix
944
 
      coarse.reinit(coarse_cell);
945
 
      FullMatrix<number> A(n, n);
946
 
      for (unsigned int k=0;k<nq;++k)
947
 
        for (unsigned int i=0;i<n;++i)
948
 
          for (unsigned int j=0;j<n;++j)
949
 
            if (fe.is_primitive())
950
 
              A(i,j) += coarse.JxW(k)
951
 
                        * coarse.shape_value(i,k)
952
 
                        * coarse.shape_value(j,k);
953
 
            else
954
 
              for (unsigned int d=0;d<nd;++d)
955
 
                A(i,j) = coarse.JxW(k)
956
 
                         * coarse.shape_value_component(i,k,d)
957
 
                         * coarse.shape_value_component(j,k,d);
958
 
  
959
 
      Householder<double> H(A);
960
1005
 
961
1006
      Vector<number> v_coarse(n);
962
1007
      Vector<number> v_fine(n);
963
 
  
 
1008
 
964
1009
      for (unsigned int cell_number=0;cell_number<nc;++cell_number)
965
1010
        {
966
1011
          FullMatrix<double> &this_matrix = matrices[ref_case-1][cell_number];
967
 
      
 
1012
 
968
1013
                                           // Compute right hand side,
969
1014
                                           // which is a fine level basis
970
1015
                                           // function tested with the
974
1019
                                    fine.get_JxW_values());
975
1020
          FEValues<dim> coarse (mapping, fe, q_coarse, update_values);
976
1021
          coarse.reinit(coarse_cell);
977
 
      
 
1022
 
978
1023
                                           // Build RHS
979
1024
 
 
1025
          const std::vector<double> & JxW = fine.get_JxW_values();
 
1026
 
980
1027
                                           // Outer loop over all fine
981
1028
                                           // grid shape functions phi_j
982
1029
          for (unsigned int j=0;j<fe.dofs_per_cell;++j)
983
1030
            {
984
 
              v_fine = 0.;
985
 
                                               // Loop over all quadrature points
986
 
              for (unsigned int k=0;k<fine.n_quadrature_points;++k)
 
1031
              for (unsigned int i=0; i<fe.dofs_per_cell;++i)
987
1032
                {
988
 
                                                   // integrate the scalar
989
 
                                                   // product
990
 
                                                   // (phi_i,phi_j) for
991
 
                                                   // all coarse shape
992
 
                                                   // functions to get the
993
 
                                                   // right hand side
994
 
                  for (unsigned int i=0;i<fe.dofs_per_cell;++i)
995
 
                    {
996
 
                      if (fe.is_primitive())
997
 
                        v_fine(i) += fine.JxW(k)
998
 
                                     * coarse.shape_value(i,k)
999
 
                                     * fine.shape_value(j,k);
1000
 
                      else
1001
 
                        for (unsigned int d=0;d<nd;++d)
1002
 
                          v_fine(i) += fine.JxW(k)
1003
 
                                       * coarse.shape_value_component(i,k,d)
1004
 
                                       * fine.shape_value_component(j,k,d);
 
1033
                  if (fe.is_primitive())
 
1034
                    {
 
1035
                      const double * coarse_i = &coarse.shape_value(i,0);
 
1036
                      const double * fine_j = &fine.shape_value(j,0);
 
1037
 
 
1038
                      double update = 0;
 
1039
                      for (unsigned int k=0; k<nq; ++k)
 
1040
                        update += JxW[k] * coarse_i[k] * fine_j[k];
 
1041
                      v_fine(i) = update;
 
1042
                    }
 
1043
                  else
 
1044
                    {
 
1045
                      double update = 0;
 
1046
                      for (unsigned int d=0; d<nd; ++d)
 
1047
                        for (unsigned int k=0; k<nq; ++k)
 
1048
                          update += JxW[k] * coarse.shape_value_component(i,k,d)
 
1049
                                           * fine.shape_value_component(j,k,d);
 
1050
                      v_fine(i) = update;
1005
1051
                    }
1006
1052
                }
 
1053
 
1007
1054
                                               // RHS ready. Solve system
1008
1055
                                               // and enter row into
1009
1056
                                               // matrix
1010
 
              H.least_squares(v_coarse, v_fine);
 
1057
              mass.vmult (v_coarse, v_fine);
1011
1058
              for (unsigned int i=0;i<fe.dofs_per_cell;++i)
1012
1059
                this_matrix(i,j) = v_coarse(i);
1013
1060
            }
1014
 
      
 
1061
 
1015
1062
                                           // Remove small entries from
1016
1063
                                           // the matrix
1017
1064
          for (unsigned int i=0; i<this_matrix.m(); ++i)
1036
1083
  ConstraintMatrix dummy;
1037
1084
  dummy.close();
1038
1085
  interpolate(dof1, u1, dof2, dummy, u2);
1039
 
}  
 
1086
}
1040
1087
 
1041
1088
 
1042
1089
 
1043
1090
template <int dim, int spacedim,
1044
1091
          template <int, int> class DH1,
1045
 
          template <int, int> class DH2,              
 
1092
          template <int, int> class DH2,
1046
1093
          class InVector, class OutVector>
1047
1094
void
1048
1095
FETools::interpolate (const DH1<dim, spacedim> &dof1,
1054
1101
  Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch());
1055
1102
 
1056
1103
  Assert(u1.size()==dof1.n_dofs(),
1057
 
         ExcDimensionMismatch(u1.size(), dof1.n_dofs()));  
 
1104
         ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
1058
1105
  Assert(u2.size()==dof2.n_dofs(),
1059
1106
         ExcDimensionMismatch(u2.size(), dof2.n_dofs()));
1060
1107
 
1073
1120
           std::map<const FiniteElement<dim,spacedim> *,
1074
1121
                    std_cxx1x::shared_ptr<FullMatrix<double> > > >
1075
1122
    interpolation_matrices;
1076
 
  
 
1123
 
1077
1124
  typename DH1<dim,spacedim>::active_cell_iterator cell1 = dof1.begin_active(),
1078
1125
                                                   endc1 = dof1.end();
1079
1126
  typename DH2<dim,spacedim>::active_cell_iterator cell2 = dof2.begin_active(),
1084
1131
  dofs.reserve (DoFTools::max_dofs_per_cell (dof2));
1085
1132
  u2 = 0;
1086
1133
 
1087
 
  for (; cell1!=endc1; ++cell1, ++cell2) 
 
1134
  for (; cell1!=endc1; ++cell1, ++cell2)
1088
1135
    {
1089
1136
      Assert(cell1->get_fe().n_components() == cell2->get_fe().n_components(),
1090
1137
             ExcDimensionMismatch (cell1->get_fe().n_components(),
1100
1147
      const bool hanging_nodes_not_allowed
1101
1148
        = ((cell2->get_fe().dofs_per_vertex != 0) &&
1102
1149
           (constraints.n_constraints() == 0));
1103
 
  
 
1150
 
1104
1151
      if (hanging_nodes_not_allowed)
1105
1152
        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1106
1153
          Assert (cell1->at_boundary(face) ||
1107
1154
                  cell1->neighbor(face)->level() == cell1->level(),
1108
1155
                  ExcHangingNodesNotAllowed(0));
1109
 
      
1110
 
      
 
1156
 
 
1157
 
1111
1158
      const unsigned int dofs_per_cell1 = cell1->get_fe().dofs_per_cell;
1112
1159
      const unsigned int dofs_per_cell2 = cell2->get_fe().dofs_per_cell;
1113
1160
      u1_local.reinit (dofs_per_cell1);
1117
1164
                                       // matrix for this particular
1118
1165
                                       // pair of elements is already
1119
1166
                                       // there
1120
 
      if (interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()] == 0)
 
1167
      if (interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()].get() ==
 
1168
          0)
1121
1169
        {
1122
1170
          std_cxx1x::shared_ptr<FullMatrix<double> >
1123
1171
            interpolation_matrix (new FullMatrix<double> (dofs_per_cell2,
1124
1172
                                                          dofs_per_cell1));
1125
1173
          interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()]
1126
1174
            = interpolation_matrix;
1127
 
            
 
1175
 
1128
1176
          FETools::get_interpolation_matrix(cell1->get_fe(),
1129
1177
                                            cell2->get_fe(),
1130
1178
                                            *interpolation_matrix);
1131
1179
        }
1132
 
      
 
1180
 
1133
1181
      cell1->get_dof_values(u1, u1_local);
1134
1182
      interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()]
1135
1183
        ->vmult(u2_local, u1_local);
1171
1219
{
1172
1220
  Assert(dof1.get_fe().n_components() == fe2.n_components(),
1173
1221
         ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components()));
1174
 
  Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));  
 
1222
  Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
1175
1223
  Assert(u1_interpolated.size()==dof1.n_dofs(),
1176
1224
         ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs()));
1177
 
  
 
1225
 
1178
1226
                                   // For continuous elements on grids
1179
1227
                                   // with hanging nodes we need
1180
1228
                                   // hanging node
1189
1237
 
1190
1238
  Vector<typename OutVector::value_type> u1_local(dofs_per_cell1);
1191
1239
  Vector<typename OutVector::value_type> u1_int_local(dofs_per_cell1);
1192
 
  
 
1240
 
1193
1241
  typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof1.begin_active(),
1194
1242
                                                 endc = dof1.end();
1195
1243
 
1196
1244
  FullMatrix<double> interpolation_matrix(dofs_per_cell1, dofs_per_cell1);
1197
1245
  FETools::get_back_interpolation_matrix(dof1.get_fe(), fe2,
1198
1246
                                         interpolation_matrix);
1199
 
  for (; cell!=endc; ++cell) 
 
1247
  for (; cell!=endc; ++cell)
1200
1248
    {
1201
1249
      if (hanging_nodes_not_allowed)
1202
1250
        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1211
1259
}
1212
1260
 
1213
1261
 
1214
 
  
 
1262
 
1215
1263
template <int dim,
1216
1264
          template <int> class DH,
1217
1265
          class InVector, class OutVector, int spacedim>
1222
1270
                          OutVector                &u1_interpolated)
1223
1271
{
1224
1272
  Assert(u1.size() == dof1.n_dofs(),
1225
 
         ExcDimensionMismatch(u1.size(), dof1.n_dofs()));  
 
1273
         ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
1226
1274
  Assert(u1_interpolated.size() == dof1.n_dofs(),
1227
1275
         ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs()));
1228
 
  
 
1276
 
1229
1277
  Vector<typename OutVector::value_type> u1_local(DoFTools::max_dofs_per_cell(dof1));
1230
1278
  Vector<typename OutVector::value_type> u1_int_local(DoFTools::max_dofs_per_cell(dof1));
1231
 
  
 
1279
 
1232
1280
  typename DH<dim>::active_cell_iterator cell = dof1.begin_active(),
1233
1281
                                         endc = dof1.end();
1234
1282
 
1237
1285
                                   // matrices
1238
1286
  std::map<const FiniteElement<dim> *,
1239
1287
           std_cxx1x::shared_ptr<FullMatrix<double> > > interpolation_matrices;
1240
 
  
1241
 
  for (; cell!=endc; ++cell) 
 
1288
 
 
1289
  for (; cell!=endc; ++cell)
1242
1290
    {
1243
1291
      Assert(cell->get_fe().n_components() == fe2.n_components(),
1244
1292
             ExcDimensionMismatch(cell->get_fe().n_components(),
1245
1293
                                  fe2.n_components()));
1246
 
      
 
1294
 
1247
1295
                                       // For continuous elements on
1248
1296
                                       // grids with hanging nodes we
1249
1297
                                       // need hanging node
1253
1301
                                       // constraints are allowed.
1254
1302
      const bool hanging_nodes_not_allowed=
1255
1303
        (cell->get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0);
1256
 
      
 
1304
 
1257
1305
      if (hanging_nodes_not_allowed)
1258
1306
        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1259
1307
          Assert (cell->at_boundary(face) ||
1261
1309
                  ExcHangingNodesNotAllowed(0));
1262
1310
 
1263
1311
      const unsigned int dofs_per_cell1 = cell->get_fe().dofs_per_cell;
1264
 
      
 
1312
 
1265
1313
                                       // make sure back_interpolation
1266
1314
                                       // matrix is available
1267
1315
      if (interpolation_matrices[&cell->get_fe()] != 0)
1271
1319
            (new FullMatrix<double>(dofs_per_cell1, dofs_per_cell1));
1272
1320
          get_back_interpolation_matrix(dof1.get_fe(), fe2,
1273
1321
                                        *interpolation_matrices[&cell->get_fe()]);
1274
 
        }      
1275
 
      
 
1322
        }
 
1323
 
1276
1324
      u1_local.reinit (dofs_per_cell1);
1277
1325
      u1_int_local.reinit (dofs_per_cell1);
1278
 
      
 
1326
 
1279
1327
      cell->get_dof_values(u1, u1_local);
1280
1328
      interpolation_matrices[&cell->get_fe()]->vmult(u1_int_local, u1_local);
1281
1329
      cell->set_dof_values(u1_int_local, u1_interpolated);
1283
1331
}
1284
1332
 
1285
1333
 
1286
 
  
 
1334
 
1287
1335
template <int dim, class InVector, class OutVector, int spacedim>
1288
1336
void FETools::back_interpolate(const DoFHandler<dim,spacedim> &dof1,
1289
1337
                               const ConstraintMatrix &constraints1,
1303
1351
    {
1304
1352
      Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(),
1305
1353
             ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components()));
1306
 
      Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));  
 
1354
      Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
1307
1355
      Assert(u1_interpolated.size()==dof1.n_dofs(),
1308
1356
             ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs()));
1309
 
      
 
1357
 
1310
1358
                                       // For continuous elements
1311
1359
                                       // first interpolate to dof2,
1312
1360
                                       // taking into account
1321
1369
}
1322
1370
 
1323
1371
 
1324
 
  
 
1372
 
1325
1373
template <int dim, class InVector, class OutVector, int spacedim>
1326
1374
void FETools::interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
1327
1375
                                        const InVector &u1,
1330
1378
{
1331
1379
  Assert(dof1.get_fe().n_components() == fe2.n_components(),
1332
1380
         ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components()));
1333
 
  Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));  
 
1381
  Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
1334
1382
  Assert(u1_difference.size()==dof1.n_dofs(),
1335
1383
         ExcDimensionMismatch(u1_difference.size(), dof1.n_dofs()));
1336
 
  
 
1384
 
1337
1385
                                   // For continuous elements on grids
1338
1386
                                   // with hanging nodes we need
1339
1387
                                   // hnaging node
1348
1396
 
1349
1397
  Vector<typename OutVector::value_type> u1_local(dofs_per_cell);
1350
1398
  Vector<typename OutVector::value_type> u1_diff_local(dofs_per_cell);
1351
 
  
 
1399
 
1352
1400
  FullMatrix<double> difference_matrix(dofs_per_cell, dofs_per_cell);
1353
1401
  FETools::get_interpolation_difference_matrix(dof1.get_fe(), fe2,
1354
1402
                                               difference_matrix);
1355
 
  
 
1403
 
1356
1404
  typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof1.begin_active(),
1357
1405
                                                 endc = dof1.end();
1358
 
  
 
1406
 
1359
1407
  for (; cell!=endc; ++cell)
1360
1408
    {
1361
1409
      if (hanging_nodes_not_allowed)
1394
1442
    }
1395
1443
}
1396
1444
 
1397
 
  
 
1445
 
1398
1446
template <int dim, class InVector, class OutVector, int spacedim>
1399
1447
void FETools::project_dg(const DoFHandler<dim,spacedim> &dof1,
1400
1448
                         const InVector &u1,
1404
1452
  Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch());
1405
1453
  Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(),
1406
1454
         ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components()));
1407
 
  Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));  
 
1455
  Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
1408
1456
  Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs()));
1409
1457
 
1410
1458
  typename DoFHandler<dim,spacedim>::active_cell_iterator cell1 = dof1.begin_active();
1413
1461
 
1414
1462
  const unsigned int n1 = dof1.get_fe().dofs_per_cell;
1415
1463
  const unsigned int n2 = dof2.get_fe().dofs_per_cell;
1416
 
  
 
1464
 
1417
1465
  Vector<double> u1_local(n1);
1418
1466
  Vector<double> u2_local(n2);
1419
1467
  std::vector<unsigned int> dofs(n2);
1420
 
  
 
1468
 
1421
1469
  FullMatrix<double> matrix(n2,n1);
1422
1470
  get_projection_matrix(dof1.get_fe(), dof2.get_fe(), matrix);
1423
 
  
 
1471
 
1424
1472
  while (cell2 != end)
1425
1473
    {
1426
1474
      cell1->get_dof_values(u1, u1_local);
1430
1478
        {
1431
1479
          u2(dofs[i])+=u2_local(i);
1432
1480
        }
1433
 
     
 
1481
 
1434
1482
      ++cell1;
1435
1483
      ++cell2;
1436
1484
    }
1437
1485
}
1438
1486
 
1439
 
  
 
1487
 
1440
1488
template <int dim, class InVector, class OutVector, int spacedim>
1441
1489
void FETools::extrapolate(const DoFHandler<dim,spacedim> &dof1,
1442
1490
                          const InVector &u1,
1480
1528
                                            endc = dof2.end(0);
1481
1529
    for (; cell!=endc; ++cell)
1482
1530
      Assert (cell->has_children(), ExcGridNotRefinedAtLeastOnce());
1483
 
  } 
 
1531
  }
1484
1532
 
1485
1533
                                   // then traverse grid bottom up
1486
1534
  for (unsigned int level=0; level<dof1.get_tria().n_levels()-1; ++level)
1540
1588
                                   // for this, acquire the lock
1541
1589
                                   // until we quit this function
1542
1590
  Threads::ThreadMutex::ScopedLock lock(fe_name_map_lock);
1543
 
  
 
1591
 
1544
1592
  Assert(fe_name_map.find(name) == fe_name_map.end(),
1545
1593
         ExcMessage("Cannot change existing element in finite element name list"));
1546
 
  
 
1594
 
1547
1595
                                   // Insert the normalized name into
1548
1596
                                   // the map
1549
1597
  fe_name_map[name] = FEFactoryPointer(factory);
1566
1614
        name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_"));
1567
1615
      const std::string name_part(name, 0, name_end);
1568
1616
      name.erase(0, name_part.size());
1569
 
  
 
1617
 
1570
1618
                                       // now things get a little more
1571
1619
                                       // complicated: FESystem. it's
1572
1620
                                       // more complicated, since we
1611
1659
                                                       // and read this
1612
1660
                                                       // multiplicity
1613
1661
                      name.erase(0,1);
1614
 
                  
 
1662
 
1615
1663
                      const std::pair<int,unsigned int> tmp
1616
1664
                        = Utilities::get_integer_at_position (name, 0);
1617
1665
                      name.erase(0, tmp.second);
1624
1672
                                                     // multiplicity is
1625
1673
                                                     // 1
1626
1674
                    base_multiplicities.push_back (1);
1627
 
              
 
1675
 
1628
1676
                                                   // so that's it for
1629
1677
                                                   // this base
1630
1678
                                                   // element. base
1637
1685
                                                   // '-'
1638
1686
                }
1639
1687
              while (name[0] == '-');
1640
 
          
 
1688
 
1641
1689
                                               // so we got to the end
1642
1690
                                               // of the '-' separated
1643
1691
                                               // list. make sure that
1651
1699
                      &&
1652
1700
                      (base_fes.size() > 0),
1653
1701
                      ExcInternalError());
1654
 
            
 
1702
 
1655
1703
                                               // ok, apparently
1656
1704
                                               // everything went ok. so
1657
1705
                                               // generate the composed
1674
1722
                                                       base_multiplicities[1]);
1675
1723
                    break;
1676
1724
                  }
1677
 
                
 
1725
 
1678
1726
                  case 3:
1679
1727
                  {
1680
1728
                    system_element = new FESystem<dim>(*base_fes[0],
1695
1743
                                               // any more
1696
1744
              for (unsigned int i=0; i<base_fes.size(); ++i)
1697
1745
                delete base_fes[i];
1698
 
            
 
1746
 
1699
1747
                                               // finally return our
1700
1748
                                               // findings
1701
1749
                                               // Add the closing ']' to
1702
1750
                                               // the length
1703
1751
              return system_element;
1704
 
            
 
1752
 
1705
1753
            }
1706
1754
          catch (...)
1707
1755
            {
1730
1778
                                           // Make sure no other thread
1731
1779
                                           // is just adding an element
1732
1780
          Threads::ThreadMutex::ScopedLock lock (fe_name_map_lock);
1733
 
        
 
1781
 
1734
1782
          AssertThrow (fe_name_map.find(name_part) != fe_name_map.end(),
1735
1783
                       FETools::ExcInvalidFEName(name));
1736
1784
                                           // Now, just the (degree)
1751
1799
              unsigned int position = name.find('(');
1752
1800
              const std::string quadrature_name(name, 0, position-1);
1753
1801
              name.erase(0,position);
1754
 
              if (quadrature_name.compare("QGaussLobatto") == 0)                
 
1802
              if (quadrature_name.compare("QGaussLobatto") == 0)
1755
1803
                {
1756
1804
                  const std::pair<int,unsigned int> tmp
1757
1805
                    = Utilities::get_integer_at_position (name, 0);
1760
1808
//return fe_name_map.find(name_part)->second->get(QGaussLobatto<1>(tmp.first));
1761
1809
                  AssertThrow (false, ExcNotImplemented());
1762
1810
                }
1763
 
              else 
 
1811
              else
1764
1812
                {
1765
1813
                  AssertThrow (false,ExcNotImplemented());
1766
1814
                }
1767
1815
            }
1768
1816
        }
1769
 
  
1770
 
    
 
1817
 
 
1818
 
1771
1819
                                       // hm, if we have come thus far, we
1772
1820
                                       // didn't know what to do with the
1773
1821
                                       // string we got. so do as the docs
1784
1832
 
1785
1833
 
1786
1834
 
1787
 
    
 
1835
 
1788
1836
 
1789
1837
template <int dim>
1790
1838
FiniteElement<dim, dim> *
1798
1846
       pos1 < name.size();
1799
1847
       pos1 = name.find('<'))
1800
1848
    {
1801
 
      
 
1849
 
1802
1850
      const unsigned int pos2 = name.find('>');
1803
1851
                                       // If there is only a single
1804
1852
                                       // character between those two,
1813
1861
        }
1814
1862
      else
1815
1863
        Assert(pos2-pos1 == 4, ExcInvalidFEName(name));
1816
 
      
 
1864
 
1817
1865
                                       // If pos1==pos2, then we are
1818
1866
                                       // probably at the end of the
1819
1867
                                       // string
1827
1875
       pos < name.size();
1828
1876
       pos = name.find("^dim"))
1829
1877
    name.erase(pos+2, 2);
1830
 
  
 
1878
 
1831
1879
                                   // Replace all occurences of "^d"
1832
1880
                                   // by using the actual dimension
1833
1881
  for (unsigned int pos = name.find("^d");
1834
1882
       pos < name.size();
1835
1883
       pos = name.find("^d"))
1836
1884
    name.at(pos+1) = '0' + dim;
1837
 
  
 
1885
 
1838
1886
  try
1839
 
    {    
 
1887
    {
1840
1888
      FiniteElement<dim,dim> *fe = internal::get_fe_from_name<dim,dim> (name);
1841
1889
 
1842
1890
                                       // Make sure the auxiliary function
1926
1974
 
1927
1975
 
1928
1976
 
 
1977
template <int dim, int spacedim>
 
1978
void
 
1979
FETools::
 
1980
compute_projection_from_face_quadrature_points_matrix (const FiniteElement<dim, spacedim> &fe,
 
1981
                                                const Quadrature<dim-1>    &lhs_quadrature,
 
1982
                                                const Quadrature<dim-1>    &rhs_quadrature,
 
1983
                                                const typename DoFHandler<dim, spacedim>::active_cell_iterator & cell,
 
1984
                                                unsigned int face,
 
1985
                                                FullMatrix<double>       &X)
 
1986
{
 
1987
  Assert (fe.n_components() == 1, ExcNotImplemented());
 
1988
  Assert (lhs_quadrature.size () > fe.degree, ExcNotGreaterThan (lhs_quadrature.size (), fe.degree));
 
1989
 
 
1990
 
 
1991
 
 
1992
                                   // build the matrices M and Q
 
1993
                                   // described in the documentation
 
1994
  FullMatrix<double> M (fe.dofs_per_cell, fe.dofs_per_cell);
 
1995
  FullMatrix<double> Q (fe.dofs_per_cell, rhs_quadrature.size());
 
1996
 
 
1997
  {
 
1998
                                // need an FEFaceValues object to evaluate shape function
 
1999
                                // values on the specified face.
 
2000
    FEFaceValues <dim> fe_face_values (fe, lhs_quadrature, update_values);
 
2001
    fe_face_values.reinit (cell, face); // setup shape_value on this face.
 
2002
 
 
2003
    for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
 
2004
      for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
 
2005
        for (unsigned int q=0; q<lhs_quadrature.size(); ++q)
 
2006
          M(i,j) += fe_face_values.shape_value (i, q) *
 
2007
                    fe_face_values.shape_value (j, q) *
 
2008
                    lhs_quadrature.weight(q);
 
2009
    for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
 
2010
    {
 
2011
      M(i,i) = (M(i,i) == 0 ? 1 : M(i,i));
 
2012
    }
 
2013
  }
 
2014
 
 
2015
  {
 
2016
    FEFaceValues <dim> fe_face_values (fe, rhs_quadrature, update_values);
 
2017
    fe_face_values.reinit (cell, face); // setup shape_value on this face.
 
2018
 
 
2019
    for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
 
2020
      for (unsigned int q=0; q<rhs_quadrature.size(); ++q)
 
2021
        Q(i,q) += fe_face_values.shape_value (i, q) *
 
2022
                  rhs_quadrature.weight(q);
 
2023
  }
 
2024
                                   // then invert M
 
2025
  FullMatrix<double> M_inverse (fe.dofs_per_cell, fe.dofs_per_cell);
 
2026
  M_inverse.invert (M);
 
2027
 
 
2028
                                   // finally compute the result
 
2029
  X.reinit (fe.dofs_per_cell, rhs_quadrature.size());
 
2030
  M_inverse.mmult (X, Q);
 
2031
 
 
2032
  Assert (X.m() == fe.dofs_per_cell, ExcInternalError());
 
2033
  Assert (X.n() == rhs_quadrature.size(), ExcInternalError());
 
2034
}
 
2035
 
 
2036
 
 
2037
 
 
2038
 
1929
2039
/*-------------- Explicit Instantiations -------------------------------*/
1930
2040
 
1931
2041
 
1943
2053
template
1944
2054
void FETools::compute_block_renumbering (
1945
2055
  const FiniteElement<deal_II_dimension>& element,
1946
 
  std::vector<unsigned int>&, std::vector<unsigned int>&_indices);
 
2056
  std::vector<unsigned int>&, std::vector<unsigned int>&_indices, bool);
1947
2057
template
1948
2058
void FETools::get_interpolation_matrix<deal_II_dimension>
1949
2059
(const FiniteElement<deal_II_dimension> &,
2006
2116
 
2007
2117
#if deal_II_dimension != 3
2008
2118
template
 
2119
void FETools::compute_block_renumbering (
 
2120
  const FiniteElement<deal_II_dimension,deal_II_dimension+1>& element,
 
2121
  std::vector<unsigned int>&, std::vector<unsigned int>&_indices, bool);
 
2122
template
2009
2123
void FETools::interpolate<deal_II_dimension,deal_II_dimension+1>
2010
2124
(const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const Vector<double> &,
2011
2125
 const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, Vector<double> &);
2360
2474
FETools::get_fe_from_name<deal_II_dimension> (const std::string &);
2361
2475
 
2362
2476
 
2363
 
template 
 
2477
template
2364
2478
void FETools::add_fe_name<deal_II_dimension>(
2365
2479
  const std::string& name,
2366
2480
  const FEFactoryBase<deal_II_dimension>* factory);
2367
 
  
 
2481
 
2368
2482
template
2369
2483
void
2370
2484
FETools::
2380
2494
                                                   const Quadrature<deal_II_dimension>    &quadrature,
2381
2495
                                                   FullMatrix<double>       &I_q);
2382
2496
 
 
2497
#if deal_II_dimension != 1
 
2498
template
 
2499
void
 
2500
FETools::
 
2501
compute_projection_from_face_quadrature_points_matrix (const FiniteElement<deal_II_dimension> &fe,
 
2502
                                                const Quadrature<deal_II_dimension-1>    &lhs_quadrature,
 
2503
                                                const Quadrature<deal_II_dimension-1>    &rhs_quadrature,
 
2504
                                                const DoFHandler<deal_II_dimension>::active_cell_iterator & cell,
 
2505
                                                unsigned int face,
 
2506
                                                FullMatrix<double>       &X);
 
2507
#endif
 
2508
 
2383
2509
 
2384
2510
 
2385
2511
/*----------------------------   fe_tools.cc     ---------------------------*/