[Getfem-commits] r5131 - in /trunk/getfem: interface/tests/matlab/ tests/

[Yves . Renard]

Author: renard
Date: Wed Nov 11 13:56:23 2015
New Revision: 5131

URL: http://svn.gna.org/viewcvs/getfem?rev=5131&view=rev
Log:
minor changes

Modified:
    trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m
    trunk/getfem/tests/laplacian.cc
    trunk/getfem/tests/laplacian.param
    trunk/getfem/tests/laplacian_with_bricks.cc
    trunk/getfem/tests/laplacian_with_bricks.param
    trunk/getfem/tests/laplacian_with_bricks.pl

Modified: trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m?rev=5131&r1=5130&r2=5131&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m     (original)
+++ trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m     Wed Nov 11 
13:56:23 2015
@@ -29,7 +29,7 @@
 r = 1e8;         % Penalization parameter for the Dirichlet condition
 draw = false;
 quadrangles = true;
-NX = 150;
+NX = 20;
 K = 2;           % Degree of the discontinuous finite element method
 interior_penalty_factor = 300*NX; % Parameter of the interior penalty term
 verify_neighbour_computation = true;
@@ -126,10 +126,6 @@
 % gf_model_set(md, 'add linear generic assembly brick', mim, 
sprintf('alpha*((%s).(%s))', jump, test_jump), INNER_FACES);
 gf_model_set(md, 'add linear generic assembly brick', mim, 
sprintf('-((%s).(%s))-((%s).(%s))+alpha*((%s).(%s))', grad_mean, test_jump, 
jump, grad_test_mean, jump, test_jump), INNER_FACES);
 
-tic;
-gf_model_get(md, 'assembly');
-toc
-
 gf_model_get(md, 'solve', 'noisy');
 U = gf_model_get(md, 'variable', 'u');
 

Modified: trunk/getfem/tests/laplacian.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/tests/laplacian.cc?rev=5131&r1=5130&r2=5131&view=diff
==============================================================================
--- trunk/getfem/tests/laplacian.cc     (original)
+++ trunk/getfem/tests/laplacian.cc     Wed Nov 11 13:56:23 2015
@@ -93,7 +93,7 @@
                              */
   std::string datafilename;
   bgeot::md_param PARAM;
-
+  
   void assembly(void);
   bool solve(void);
   void init(void);
@@ -111,6 +111,7 @@
   std::string FEM_TYPE  = PARAM.string_value("FEM_TYPE","FEM name");
   std::string INTEGRATION = PARAM.string_value("INTEGRATION",
                                               "Name of integration method");
+
   cout << "MESH_TYPE=" << MESH_TYPE << "\n";
   cout << "FEM_TYPE="  << FEM_TYPE << "\n";
   cout << "INTEGRATION=" << INTEGRATION << "\n";

Modified: trunk/getfem/tests/laplacian.param
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/tests/laplacian.param?rev=5131&r1=5130&r2=5131&view=diff
==============================================================================
--- trunk/getfem/tests/laplacian.param  (original)
+++ trunk/getfem/tests/laplacian.param  Wed Nov 11 13:56:23 2015
@@ -8,7 +8,7 @@
 FT = 0.1;                   % parameter for the exact solution.
 N = 2;                      % domain in dimension N
 QUAD = 0;                   % use quadrilaterons or not
-NX = 100;                    % space step.
+NX = 100;                   % space step.
 MESH_NOISED = 0;            % Set to one if you want to "shake" the mesh
 
 if (N == 1)

Modified: trunk/getfem/tests/laplacian_with_bricks.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/tests/laplacian_with_bricks.cc?rev=5131&r1=5130&r2=5131&view=diff
==============================================================================
--- trunk/getfem/tests/laplacian_with_bricks.cc (original)
+++ trunk/getfem/tests/laplacian_with_bricks.cc Wed Nov 11 13:56:23 2015
@@ -73,7 +73,7 @@
 */
 struct laplacian_problem {
 
-  enum { DIRICHLET_BOUNDARY_NUM = 0, NEUMANN_BOUNDARY_NUM = 1};
+  enum { DIRICHLET_BOUNDARY_NUM = 0, NEUMANN_BOUNDARY_NUM = 1, INNER_FACES = 
2};
   enum { DIRICHLET_WITH_MULTIPLIERS = 0, DIRICHLET_WITH_PENALIZATION = 1};
   getfem::mesh mesh;        /* the mesh */
   getfem::mesh_im mim;      /* the integration methods. */
@@ -87,6 +87,8 @@
 
   std::string datafilename;
   bgeot::md_param PARAM;
+  bool DG_TERMS;
+  scalar_type interior_penalty_factor;
 
   bool solve(void);
   void init(void);
@@ -103,6 +105,11 @@
   std::string FEM_TYPE  = PARAM.string_value("FEM_TYPE","FEM name");
   std::string INTEGRATION = PARAM.string_value("INTEGRATION",
                                               "Name of integration method");
+  DG_TERMS = (PARAM.int_value("DG_TERMS",
+                             "Discontinuous Galerkin terms or not") != 0);
+  interior_penalty_factor = PARAM.real_value("INTERIOR_PENALTY_FACTOR",
+                                            "Interior penalty factor for DG");
+
   cout << "MESH_TYPE=" << MESH_TYPE << "\n";
   cout << "FEM_TYPE="  << FEM_TYPE << "\n";
   cout << "INTEGRATION=" << INTEGRATION << "\n";
@@ -176,6 +183,10 @@
       mesh.region(DIRICHLET_BOUNDARY_NUM).add(i.cv(), i.f());
     }
   }
+
+  if (DG_TERMS)
+    mesh.region(INNER_FACES) = getfem::inner_faces_of_mesh(mesh);
+
 }
 
 bool laplacian_problem::solve(void) {
@@ -215,8 +226,22 @@
       (model, mim, "u", dirichlet_coefficient,
        DIRICHLET_BOUNDARY_NUM, "DirichletData");
   
+  if (DG_TERMS) {
+    model.add_initialized_scalar_data("alpha", interior_penalty_factor);
+    std::string jump="((u-Interpolate(u,neighbour_elt))*Normal)";
+    std::string 
test_jump="((Test_u-Interpolate(Test_u,neighbour_elt))*Normal)";
+    std::string grad_mean="((Grad_u+Interpolate(Grad_u,neighbour_elt))*0.5)";
+    std::string grad_test_mean
+      ="((Grad_Test_u+Interpolate(Grad_Test_u,neighbour_elt))*0.5)";
+    std::string expr = "-("+grad_mean+").("+test_jump+") "
+      "- ("+jump+").("+grad_test_mean+")"
+      "+ alpha*("+jump+").("+test_jump+")";
+    getfem::add_linear_generic_assembly_brick(model, mim, expr,
+                                             INNER_FACES, true);
+  }
+
   model.listvar(cout);
- 
+
   gmm::iteration iter(residual, 1, 40000);
   getfem::standard_solve(model, iter);
 

Modified: trunk/getfem/tests/laplacian_with_bricks.param
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/tests/laplacian_with_bricks.param?rev=5131&r1=5130&r2=5131&view=diff
==============================================================================
--- trunk/getfem/tests/laplacian_with_bricks.param      (original)
+++ trunk/getfem/tests/laplacian_with_bricks.param      Wed Nov 11 13:56:23 2015
@@ -8,16 +8,23 @@
 FT = 0.1;                           % parameter for the exact solution.
 N = 2;                              % domain in dimension N
 QUAD = 0;                           % use quadrilaterons or not
-NX = 100;                           % space step.
+NX = 50;                            % space step.
 MESH_NOISED = 0;                    % "shake" the mesh or not
+DG_TERMS = 0;                        % Discontinuous Galerkin terms or not
+INTERIOR_PENALTY_FACTOR = 1000;      % Interior penalty factor for DG
 DIRICHLET_VERSION = 1;              % 0 = With Lagrange multipliers
                                     % 1 = penalization.
 DIRICHLET_COEFFICIENT = 1E10;       % Penalization coefficient.
 
+
 if (N == 1)
   FT = 5.0;
   MESH_TYPE = 'GT_PK(1,1)';          % segments
-  FEM_TYPE = 'FEM_PK(1,2)';          % P1 for segments
+  if (DG_TERMS)
+    FEM_TYPE = 'FEM_PK_DISCONTINUOUS(1,2)';          % P2 for segments
+  else
+    FEM_TYPE = 'FEM_PK(1,2)';          % P2 for segments
+  end
 %  FEM_TYPE = 'FEM_HERMITE(1)'; DATA_FEM_TYPE = 'FEM_PK(1,3)';
   INTEGRATION = 'IM_GAUSS1D(6)';     % Gauss-Legendre integration on segments
 end
@@ -25,7 +32,11 @@
 if (N == 2 && ~QUAD)
   FT = 10.0;
   MESH_TYPE = 'GT_PK(2,1)';          % triangles
-  FEM_TYPE = 'FEM_PK(2,1)';          % P1 for triangles
+  if (DG_TERMS)
+    FEM_TYPE = 'FEM_PK_DISCONTINUOUS(2,2)';        % P2 for triangles
+  else
+    FEM_TYPE = 'FEM_PK(2,2)';          % P2 for triangles
+  end
   %FEM_TYPE = 'FEM_REDUCED_HCT_TRIANGLE'; DATA_FEM_TYPE = 'FEM_PK(2,3)';
   %INTEGRATION = 'IM_HCT_COMPOSITE(IM_TRIANGLE(6))';
   %FEM_TYPE = 'FEM_HERMITE(2)'; DATA_FEM_TYPE = 'FEM_PK(2,2)';
@@ -41,13 +52,21 @@
   else
     MESH_TYPE = 'GT_LINEAR_QK(2)';   % quadrilaterons (linear transformation)
   end
-  FEM_TYPE = 'FEM_QK(2,2)';          % Q1 for quadrilaterons
+  if (DG_TERMS)
+    FEM_TYPE = 'FEM_QK_DISCONTINUOUS(2,2)';          % Q2 for quadrilaterons
+  else
+    FEM_TYPE = 'FEM_QK(2,2)';          % Q2 for quadrilaterons
+  end
   INTEGRATION = 'IM_GAUSS_PARALLELEPIPED(2,6)';
 end
 
 if (N == 3 && ~QUAD)
   MESH_TYPE = 'GT_PK(3,1)';          % tetrahedrons
-  FEM_TYPE = 'FEM_PK(3,2)';          % P1 for tetrahedrons
+  if (DG_TERMS)
+    FEM_TYPE = 'FEM_PK_DISCONTINUOUS(3,2)';          % P1 for tetrahedrons
+  else
+    FEM_TYPE = 'FEM_PK(3,2)';          % P1 for tetrahedrons
+  end
   INTEGRATION = 'IM_TETRAHEDRON(2)'; % quadrature rule for polynomials up
                                      % to degree 6 on tetrahedrons
 end
@@ -58,7 +77,11 @@
   else
     MESH_TYPE = 'GT_LINEAR_QK(3)';   % hexahedrons (linear transformation)
   end
-  FEM_TYPE = 'FEM_QK(3,1)';          % Q1 for hexahedrons
+  if (DG_TERMS)
+    FEM_TYPE = 'FEM_QK_DISCONTINUOUS(3,2)';          % Q2 for hexahedrons
+  else
+    FEM_TYPE = 'FEM_QK(3,2)';          % Q2 for hexahedrons
+  end
   INTEGRATION = 'IM_GAUSS_PARALLELEPIPED(3,6)';
 end
 

Modified: trunk/getfem/tests/laplacian_with_bricks.pl
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/tests/laplacian_with_bricks.pl?rev=5131&r1=5130&r2=5131&view=diff
==============================================================================
--- trunk/getfem/tests/laplacian_with_bricks.pl (original)
+++ trunk/getfem/tests/laplacian_with_bricks.pl Wed Nov 11 13:56:23 2015
@@ -39,6 +39,8 @@
 DIRICHLET_VERSION = 1;           % 0 = With Lagrange multipliers
                                  % 1 = penalization.
 DIRICHLET_COEFFICIENT = 1E10;    % Penalization coefficient.
+DG_TERMS = 0;                     % No discontinous Galerkin terms
+INTERIOR_PENALTY_FACTOR = 0;      % Interior penalty factor for DG
 
 ;
 close(TMPF);