[Getfem-commits] r4902 - in /trunk/getfem/interface: src/gf_mesh_get.cc tests/python/demo_navier_stokes.py

[logari81]

Author: logari81
Date: Wed Mar 25 12:17:29 2015
New Revision: 4902

URL: http://svn.gna.org/viewcvs/getfem?rev=4902&view=rev
Log:
add new Mesh.outer_faces_in_box() method and incompressible Navier-Stokes demo 
(from fenics)

Added:
    trunk/getfem/interface/tests/python/demo_navier_stokes.py
Modified:
    trunk/getfem/interface/src/gf_mesh_get.cc

Modified: trunk/getfem/interface/src/gf_mesh_get.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_mesh_get.cc?rev=4902&r1=4901&r2=4902&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_mesh_get.cc   (original)
+++ trunk/getfem/interface/src/gf_mesh_get.cc   Wed Mar 25 12:17:29 2015
@@ -157,19 +157,30 @@
 
 
 static void
-outer_faces(const getfem::mesh &m, mexargs_in &in, mexargs_out &out, bool 
with_normal)
+outer_faces(const getfem::mesh &m, mexargs_in &in, mexargs_out &out, const 
std::string &condition="")
 {
   mesh_faces_by_pts_list lst;
-  dal::bit_vector convex_tested;
-  dal::bit_vector cvlst;
+  dal::bit_vector cvlst, checked_pids, rejected_pids;
+
+  bool with_normal(condition == "direction");
+  bool with_box(condition == "box");
 
   darray normal_vector;
-  bgeot::base_node un;
+  bgeot::base_node un, pmin, pmax;
   scalar_type threshold(0);
   if (with_normal) {
     normal_vector = in.pop().to_darray();
     scalar_type angle = in.pop().to_scalar();
     threshold = gmm::vect_norm2(normal_vector)*cos(angle);
+  } else if (with_box) {
+    darray p1 = in.pop().to_darray(m.dim());
+    darray p2 = in.pop().to_darray(m.dim());
+    pmin.resize(m.dim());
+    pmax.resize(m.dim());
+    for (size_type k=0; k < m.dim(); ++k) {
+      pmin[k] = std::min(p1[k],p2[k]);
+      pmax[k] = std::max(p1[k],p2[k]);
+    }
   }
 
   if (in.remaining()) cvlst = in.pop().to_bit_vector(&m.convex_index());
@@ -196,15 +207,30 @@
   size_type fcnt = 0;
   for (size_type i = 0; i < lst.size(); i++) {
     if (lst[i].cnt == 1) {
-      size_type icv = lst[i].cv;
-      short_type iif = lst[i].f;
-      if (with_normal) un = m.mean_normal_of_face_of_convex(icv, iif);
-      if (!with_normal ||
-          (gmm::vect_sp(normal_vector, un) - threshold >= -1E-8)) {
-        fcnt++;
-      } else {
-        lst[i].cnt = 2;
+      if (with_normal) {
+        un = m.mean_normal_of_face_of_convex(lst[i].cv, lst[i].f);
+        if (gmm::vect_sp(normal_vector, un) < threshold - 1E-8) {
+          lst[i].cnt = -1;
+        }
+      } else if (with_box) {
+        for (std::vector<size_type>::const_iterator pid=lst[i].ptid.begin();
+             pid != lst[i].ptid.end(); ++pid) {
+          if (!checked_pids.is_in(*pid)) {
+            checked_pids.add(*pid);
+            for (size_type k=0; k < m.dim(); ++k)
+              if (m.points()[*pid][k] < pmin[k] ||
+                  m.points()[*pid][k] > pmax[k]) {
+                rejected_pids.add(*pid);
+                break;
+              }
+          }
+          if (rejected_pids.is_in(*pid)) {
+            lst[i].cnt = -1;
+            break;
+          }
+        }
       }
+      if (lst[i].cnt == 1) fcnt++;
     }
 
   }
@@ -723,33 +749,45 @@
     /address@hidden CVFIDs = ('outer faces'[, CVIDs])
     Return the set of faces not shared by two convexes.
 
-    `CVFIDs` is a two-rows matrix, the first row lists convex #ids,
-    and the second lists face numbers (local number in the convex).
-    If `CVIDs` is not given, all convexes are considered, and it
-    basically returns the mesh boundary. If `CVIDs` is given, it
-    returns the boundary of the convex set whose #ids are listed
-    in address@hidden/
+    The output `CVFIDs` is a two-rows matrix, the first row lists
+    convex #ids, and the second one lists face numbers (local number
+    in the convex). If `CVIDs` is not given, all convexes are
+    considered, and it basically returns the mesh boundary. If `CVIDs`
+    is given, it returns the boundary of the convex set whose #ids are
+    listed in address@hidden/
     sub_command
       ("outer faces", 0, 1, 0, 1,
        check_empty_mesh(pmesh);
-       outer_faces(*pmesh, in, out, false);
+       outer_faces(*pmesh, in, out);
        );
 
     /address@hidden CVFIDs = ('outer faces with direction', @vec v, @scalar 
angle [, CVIDs])
-    Return the set of faces not shared by two convexes and whose mean outward 
unit vector have at most an angle `angle` with vector `v`.
-
-    `CVFIDs` is a two-rows matrix, the first row lists convex #ids,
-    and the second lists face numbers (local number in the convex).
-    If `CVIDs` is not given, all convexes are considered, and it
-    basically returns the mesh boundary. If `CVIDs` is given, it
-    returns the boundary of the convex set whose #ids are listed
-    in address@hidden/
+    Return the set of faces not shared by two convexes and with a mean outward 
vector lying within an angle `angle` (in radians) from vector `v`.
+
+    The output `CVFIDs` is a two-rows matrix, the first row lists convex
+    #ids, and the second one lists face numbers (local number in the
+    convex). If `CVIDs` is given, it returns portion of the boundary of
+    the convex set defined by the #ids listed in address@hidden/
     sub_command
       ("outer faces with direction", 2, 3, 0, 1,
        check_empty_mesh(pmesh);
-       outer_faces(*pmesh, in, out, true);
-       );
+       outer_faces(*pmesh, in, out, "direction");
+       );
+
     
+    /address@hidden CVFIDs = ('outer faces in box', @vec pmin, @vec pmax [, 
CVIDs])
+    Return the set of faces not shared by two convexes and lying within the 
box defined by the corner points `pmin` and `pmax`.
+
+    The output `CVFIDs` is a two-rows matrix, the first row lists convex
+    #ids, and the second one lists face numbers (local number in the
+    convex). If `CVIDs` is given, it returns portion of the boundary of
+    the convex set defined by the #ids listed in address@hidden/
+    sub_command
+      ("outer faces in box", 2, 3, 0, 1,
+       check_empty_mesh(pmesh);
+       outer_faces(*pmesh, in, out, "box");
+       );
+
     
     /address@hidden CVFIDs = ('faces from cvid'[, @ivec CVIDs][, 'merge'])
     Return a list of convexes faces from a list of convex #id.

Added: trunk/getfem/interface/tests/python/demo_navier_stokes.py
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/python/demo_navier_stokes.py?rev=4902&view=auto
==============================================================================
--- trunk/getfem/interface/tests/python/demo_navier_stokes.py   (added)
+++ trunk/getfem/interface/tests/python/demo_navier_stokes.py   Wed Mar 25 
12:17:29 2015
@@ -0,0 +1,138 @@
+#!/usr/bin/env python
+# -*- coding: UTF8 -*-
+# Python GetFEM++ interface
+#
+# Copyright (C) 2015 Konstantinos Poulios.
+#
+# This file is a part of GetFEM++
+#
+# Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+# under  the  terms  of the  GNU  Lesser General Public License as published
+# by  the  Free Software Foundation;  either version 3 of the License,  or
+# (at your option) any later version along with the GCC Runtime Library
+# Exception either version 3.1 or (at your option) any later version.
+# This program  is  distributed  in  the  hope  that it will be useful,  but
+# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+# License and GCC Runtime Library Exception for more details.
+# You  should  have received a copy of the GNU Lesser General Public License
+# along  with  this program;  if not, write to the Free Software Foundation,
+# Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+#
+############################################################################
+
+""" Incompressible Navier-Stokes equation solved in an L-shaped domain.
+
+  This program is used to check that python-getfem is working. This is also
+  a good example of use of GetFEM++.
+"""
+
+# This example is ported from:
+# 
http://fenicsproject.org/documentation/dolfin/1.0.1/python/demo/pde/navier-stokes/python/documentation.html
+
+import getfem as gf
+import numpy as np
+
+
+W1 = 0.5         #       _____
+H1 = 0.5         #      |  W1 |
+W2 = 0.5         #      |     |H1
+H2 = 0.5         #      |     |_____
+NW1 = 10         # H1+H2|        W2 |
+NH1 = 10         #      |           |H2
+NW2 = 10         #      |___________|
+NH2 = 10         #          W1+W2
+
+dt = 0.01
+T = 3.
+nu = 0.01
+
+p_in_str = "np.sin(3*{0})"
+
+#Mesh and MeshRegion
+IN_RG = 1
+OUT_RG = 2
+INOUT_RG = 3
+WALL_RG = 4
+
+geotrans = "GT_QK(2,2)"
+m = gf.Mesh("import", "structured",
+            "GT='%s';ORG=[%f,%f];SIZES=[%f,%f];NSUBDIV=[%i,%i]"
+            % (geotrans, 0, H2, W1, H1, NW1, NH1))
+m.merge(gf.Mesh("import", "structured",
+                "GT='%s';ORG=[%f,%f];SIZES=[%f,%f];NSUBDIV=[%i,%i]"
+                % (geotrans, 0, 0, W1, H2, NW1, NH2)))
+m.merge(gf.Mesh("import", "structured",
+                "GT='%s';ORG=[%f,%f];SIZES=[%f,%f];NSUBDIV=[%i,%i]"
+                % (geotrans, W1, 0, W2, H2, NW2, NH2)))
+m.optimize_structure()
+
+l_rg = m.outer_faces_with_direction([-1,0], np.pi/180)
+b_rg = m.outer_faces_with_direction([0,-1], np.pi/180)
+r_rg = m.outer_faces_with_direction([1,0], np.pi/180)
+t_rg = m.outer_faces_with_direction([0,1], np.pi/180)
+in_rg = m.outer_faces_in_box([-1e-6,H1+H2-1e-6],[W1+1e-6,H1+H2+1e-6])
+out_rg = m.outer_faces_in_box([W1+W2-1e-6,-1e-6],[W1+W2+1e-6,H2+1e-6])
+m.set_region(IN_RG, in_rg)
+m.set_region(OUT_RG, out_rg)
+m.extend_region(INOUT_RG, in_rg)
+m.extend_region(INOUT_RG, out_rg)
+
+m.extend_region(WALL_RG, l_rg)
+m.extend_region(WALL_RG, b_rg)
+m.extend_region(WALL_RG, r_rg)
+m.extend_region(WALL_RG, t_rg)
+m.region_subtract(WALL_RG, INOUT_RG)
+
+#MeshFem
+mfv_ = gf.MeshFem(m, 2)
+mfv_.set_classical_fem(2)
+kept_dofs = np.setdiff1d(np.arange(mfv_.nb_basic_dof()),
+                         mfv_.basic_dof_on_region(WALL_RG))
+mfv = gf.MeshFem("partial", mfv_, kept_dofs)
+
+mfp_ = gf.MeshFem(m, 1)
+mfp_.set_classical_fem(1)
+kept_dofs = np.setdiff1d(np.arange(mfp_.nb_basic_dof()),
+                         mfp_.basic_dof_on_region(OUT_RG))
+mfp = gf.MeshFem("partial", mfp_, kept_dofs)
+
+mim = gf.MeshIm(m, 5) # 9 gauss points per quad
+
+md = gf.Model("real")
+md.add_fem_variable("v", mfv)
+md.add_fem_data("v0", mfv)
+md.add_fem_variable("p", mfp)
+md.add_fem_data("p_in", mfp_)
+md.add_initialized_data("f", [0,0])
+md.add_initialized_data("dt", [dt])
+md.add_initialized_data("nu", [nu])
+
+md.add_Dirichlet_condition_with_multipliers(mim, "p", mfp, IN_RG, "p_in")
+md.add_nonlinear_generic_assembly_brick\
+(mim, "1/dt*(v-v0).Test_v + (Grad_v0*v0).Test_v + nu*Grad_v:Grad_Test_v - 
f.Test_v")
+md.add_nonlinear_generic_assembly_brick\
+(mim, "Grad_p.Grad_Test_p + 1/dt*Trace(Grad_v)*Test_p")
+
+mmat_v = gf.asm_mass_matrix(mim, mfv)
+#mmat_v = gf.asm_generic(mim, 2, "Test2_v.Test_v", -1, "v", 1, mfv, 
np.zeros(mfv.nbdof()))
+IV = md.interval_of_variable("v")
+IV = range(IV[0],IV[0]+IV[1])
+
+t = 0
+step = 0
+while t < T+1e-8:
+   print("Solving step at t=%f" % t)
+   md.set_variable\
+   ("p_in", mfp_.eval(p_in_str.format(t), globals(), locals()).flatten("F"))
+   md.set_variable("v0", md.variable("v"))
+   md.solve("noisy", "lsolver", "mumps", "max_res", 1e-8)
+   vv = gf.asm_generic(mim, 1, "(v-dt*Grad_p).Test_v", -1, md)[IV]
+   md.set_variable("v", gf.linsolve_mumps(mmat_v, vv))
+
+   mfv.export_to_vtk("results_%i.vtk" % step,
+                     mfv, md.variable("v"), "Velocity",
+                     mfp, md.variable("p"), "Pressure")
+   t += dt
+   step += 1
+