Commit 530ab151 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

test reload

parent d4bbdc9d
......@@ -15,7 +15,7 @@ class StrongBoundary(Structure):
class fluid_problem :
def __init__(self, mesh_file_name, g, mu, rho, epsilon, strong_boundaries, notComputeEpsilon):
def __init__(self, mesh_file_name, g, mu, rho, epsilon, strong_boundaries):
nsb = len(strong_boundaries)
class Bnd :
def __init__(self, b) :
......@@ -31,7 +31,7 @@ class fluid_problem :
asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],BNDCB(Bnd(i[2]).apply)) for i in strong_boundaries])
self.asb = asb
lib.fluid_problem_new.restype = c_void_p
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), c_double(mu), c_double(rho), c_double(epsilon), nsb, asb,c_int(notComputeEpsilon)))
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), c_double(mu), c_double(rho), c_double(epsilon), nsb, asb))
if self._fp == None :
raise NameError("cannot create fluid problem")
......@@ -39,8 +39,8 @@ class fluid_problem :
if(self._fp != None) :
lib.fluid_problem_free(self._fp)
def adapt_mesh(self, lcmax, lcmin, n_el, alpha, notComputeEpsilon) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(alpha), c_int(notComputeEpsilon))
def adapt_mesh(self, lcmax, lcmin, n_el) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el))
def export(self, output_dir, t, it) :
lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))
......
......@@ -16,7 +16,7 @@ class StrongBoundary(Structure):
class fluid_problem :
def __init__(self, mesh_file_name, g, mu, rho, alpha, strong_boundaries, notComputeEpsilon):
def __init__(self, mesh_file_name, g, mu, rho, alpha, strong_boundaries):
nsb = len(strong_boundaries)
class Bnd :
def __init__(self, b) :
......@@ -32,7 +32,7 @@ class fluid_problem :
asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],BNDCB(Bnd(i[2]).apply)) for i in strong_boundaries])
self.asb = asb
lib.fluid_problem_new.restype = c_void_p
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), c_double(mu), c_double(rho), c_double(alpha), nsb, asb, c_int(notComputeEpsilon)))
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), c_double(mu), c_double(rho), c_double(alpha), nsb, asb))
if self._fp == None :
raise NameError("cannot create fluid problem")
......@@ -40,8 +40,8 @@ class fluid_problem :
if(self._fp != None) :
lib.fluid_problem_free(self._fp)
def adapt_mesh(self, lcmax, lcmin, n_el, alpha, notComputeEpsilon) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(alpha), c_int(notComputeEpsilon))
def adapt_mesh(self, lcmax, lcmin, n_el) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el))
def export(self, output_dir, t, it) :
lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))
......
......@@ -210,7 +210,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
const double *old_porosity = problem->old_porosity;
const double mu = problem->mu;
const double rho = problem->rho;
const double *epsilon = problem->epsilon;
const double epsilon = problem->epsilon;
size_t local_size = N_SF*n_fields;
size_t N_E = mesh->n_elements;
double *all_local_vector = malloc(N_E*local_size*sizeof(double));
......@@ -280,7 +280,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
dphidp += dphii[k]*dp[k];
divu += du[k][k];
}
local_vector[iphi+N_SF*DIMENSION] += jw*(epsilon[el[iphi]]*dphidp+phii*(divu+dcdt));
local_vector[iphi+N_SF*DIMENSION] += jw*(epsilon*dphidp+phii*(divu+dcdt));
for (int jphi = 0; jphi < N_SF; ++jphi) {
double phij = phi[jphi];
const double *dphij = dphi[jphi];
......@@ -312,7 +312,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
//LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*(phii*phij/dt-phij*dphii[j]*u[j]/c);
LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt-phij/c*dcdt+udtau/c);
}
LOCAL_MATRIX(P,P) += jw*epsilon[el[iphi]]*dphiidphij;
LOCAL_MATRIX(P,P) += jw*epsilon*dphiidphij;
}
}
}
......@@ -453,7 +453,7 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
}
}
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double alpha, int n_strong_boundaries, StrongBoundary *strong_boundaries, int notComputeEpsilon)
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
{
static int initialized = 0;
if (!initialized) {
......@@ -473,6 +473,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
problem->mu = mu;
problem->g = g;
problem->mesh = mesh;
problem->epsilon = epsilon;
problem->strong_boundaries = malloc(sizeof(StrongBoundary)*n_strong_boundaries);
problem->n_strong_boundaries = n_strong_boundaries;
for (int i = 0; i < n_strong_boundaries; ++i) {
......@@ -481,7 +482,6 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
problem->strong_boundaries[i].field = strong_boundaries[i].field;
}
problem->porosity = malloc(mesh->n_nodes*sizeof(double));
problem->epsilon = malloc(mesh->n_nodes*sizeof(double));
problem->old_porosity = malloc(mesh->n_nodes*sizeof(double));
problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
for (int i = 0; i < problem->mesh->n_nodes; ++i)
......@@ -495,20 +495,6 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
}
problem->node_volume = malloc(0);
fluid_problem_compute_node_volume(problem);
#if DIMENSION == 2
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha*sqrt(problem->node_volume[i])*sqrt(problem->node_volume[i])/(problem->mu/problem->rho);
}
#else
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha*cbrt(problem->node_volume[i])*cbrt(problem->node_volume[i])/(problem->mu/problem->rho);
}
#endif
if(notComputeEpsilon){
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha;
}
}
// end to remove
problem->linear_system = NULL;
#ifdef HXT_HAVE_PETSC
......@@ -531,7 +517,6 @@ void fluid_problem_free(FluidProblem *problem)
{
free(problem->solution);
free(problem->porosity);
free(problem->epsilon);
free(problem->old_porosity);
hxtLinearSystemDelete(&problem->linear_system);
free(problem->particle_uvw);
......@@ -721,7 +706,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double alpha, int notComputeEpsilon)
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
{
struct timespec tic, toc;
fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el);
......@@ -767,8 +752,6 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
free(newx);
free(problem->solution);
problem->solution = new_solution;
free(problem->epsilon);
problem->epsilon = malloc(new_mesh->n_nodes*sizeof(double));
free(problem->porosity);
problem->porosity = malloc(sizeof(double)*new_mesh->n_nodes);
free(problem->old_porosity);
......@@ -786,20 +769,6 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
problem->particle_element_id[i] = -1;
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_node_volume(problem);
#if DIMENSION == 2
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha*sqrt(problem->node_volume[i])*sqrt(problem->node_volume[i])/(problem->mu/problem->rho);
}
#else
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha*cbrt(problem->node_volume[i])*cbrt(problem->node_volume[i])/(problem->mu/problem->rho);
}
#endif
if(notComputeEpsilon){
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha;
}
}
fluid_problem_compute_porosity(problem);
hxtLinearSystemDelete(&problem->linear_system);
#ifdef HXT_HAVE_PETSC
......@@ -809,6 +778,23 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
#endif
}
void fluid_problem_after_import(FluidProblem *problem)
{
mesh_tree_free(problem->mesh_tree);
problem->mesh_tree = mesh_tree_create(problem->mesh);
for (int i = 0; i < problem->n_particles; ++i)
problem->particle_element_id[i] = -1;
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_node_volume(problem);
fluid_problem_compute_porosity(problem);
hxtLinearSystemDelete(&problem->linear_system);
#ifdef HXT_HAVE_PETSC
hxtLinearSystemCreatePETSc(&problem->linear_system,problem->mesh->n_elements,N_N,n_fields,problem->mesh->elements,"fluid");
#else
hxtLinearSystemCreateLU(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements);
#endif
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid)
{
......
......@@ -5,6 +5,12 @@
#include "mesh_find.h"
#include "hxt_linear_system.h"
#if DIMENSION==2
#define N_N 3
#else
#define N_N 4
#endif
typedef struct {
char *tag;
int field;
......@@ -12,7 +18,7 @@ typedef struct {
}StrongBoundary;
typedef struct {
double *epsilon;
double epsilon;
double rho;
double mu;
double alpha;
......@@ -46,8 +52,9 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
int fluid_problem_implicit_euler(FluidProblem *problem, double dt);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double alpha, int n_strong_boundaries, StrongBoundary *strong_boundaries, int notComputeEpsilon);
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double alpha, int notComputeEpsilon);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el);
int *fluid_problem_particle_element_id(FluidProblem *problem);
int fluid_problem_n_particles(FluidProblem *problem);
void fluid_problem_after_import(FluidProblem *problem);
#endif
......@@ -143,7 +143,7 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
mb_xml_end_element(r,&gridelement);
mb_xml_end_element(r,&fileelement);
mb_xml_reader_destroy(r);
fluid_problem_after_import(problem);
return 0;
}
......
L = 0.05;
H = 0.6;
P = 0.02;
y = 0;
lc = 0.01;
Point(1) = {-L, H, -P, lc};
Point(2) = {-L, -H, -P, lc};
Point(3) = {L, -H, -P, lc};
Point(4) = {L, H, -P, lc};
Point(5) = {-L, H, P, lc};
Point(6) = {-L, -H, P, lc};
Point(7) = {L, -H, P, lc};
Point(8) = {L, H, P, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line(9) = {5, 1};
Line(10) = {4, 8};
Line(11) = {2, 6};
Line(12) = {3, 7};
Line Loop(1) = {1:4};
Line Loop(2) = {5:8};
Line Loop(3) = {9,-4,10,8};
Line Loop(4) = {2,12,-6,-11};
Line Loop(5) = {12,7,-10,-3};
Line Loop(6) = {1,11,-5,9};
Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Plane Surface(4) = {4};
Plane Surface(5) = {5};
Plane Surface(6) = {6};
Surface Loop(1) = {1,-3,-2,4,1,-5,-2,6};
Volume(1) = {1};
Physical Surface("Z") = {1,2};
Physical Surface("X") = {5,6};
Physical Surface("Top") = {3};
Physical Surface("Bottom") = {4};
Physical Volume("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.Algorithm = 5;
Merge "lc.pos";
Field[1] = PostView;
Field[1].IView = 0;
Background Field = 1;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 0;
#!/usr/bin/env python
from marblesbag import fluid3 as fluid
from marblesbag import scontact3
import numpy as np
import os
import time
import shutil
import random
def genInitialPosition(filename, r, rout, rhop, compacity) :
p = scontact3.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","X","Z"])
N = compacity*rout**3/(r**3)
e = 2*rout/(6*N/np.pi)**(1./3.)
for x in np.arange(rout, -rout, -e):
for y in np.arange(rout, -rout, -e):
for z in np.arange(rout, -rout, -e):
R2 = x**2 + y**2 + z**2
if R2<(rout-r)**2:
p.add_particle((x+random.uniform(-e/2+r,e/2-r), y+random.uniform(-e/2+r,e/2-r), z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop);
p.position()[:, 1] += 0.52
print(len(p.position()[:, 1] ))
p.write_vtk(filename,0,0)
outputdir = "MetzgerB1"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
initialize = 0
r=154e-6
R = 3.3e-3
compacity = .20
drho = 35.4
p = scontact3.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = -9.81
rhop = 2450
nu = 1.17/1030.
rho = 1030#rhop-drho/compacity
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 1000
#numerical parameters
lcmin = 0.0004 # approx r*100 but should match the mesh size
dt = 5e-2
alpha = 1e-3
epsilon = alpha*lcmin**2 /nu
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
mu = nu*rho
genInitialPosition(outputdir, r, R, rhop, compacity)
p = scontact3.ParticleProblem()
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 10
outf1 = 100000
strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
ii = 0
jj = 0
tEnd1 = 100
print(p.velocity()[1,1])
print(1-compacity)
if initialize:
for jj in range(4):
dt1 = dt
t = 0
if jj!=0:
fluid.adapt_mesh(5e-3,8e-4,600000)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd1 :
niter = fluid.implicit_euler(dt1)
forces = fluid.compute_node_force(dt1)
vn = p.velocity() + forces * dt1 / p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
if niter < 5 and dt1 < 20:
dt1 *= 3
elif niter > 5 and dt1 > dt:
dt1 /= 3
print(dt1)
t += dt1
ii += 1
ii = 0
t = 0
tic = time.clock()
fluid.export_vtk(outputdir,0,0)
forces = g*p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
if (ii%5==0 and ii != 0):
fluid.adapt_mesh(5e-3,8e-4,600000)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
for i in range(nsub) :
p.iterate(dt/nsub, forces)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
#!/usr/bin/env python
from marblesbag import fluid3 as fluid
from marblesbag import scontact3
import numpy as np
import os
import time
import shutil
import random
outputdir = "MetzgerB1"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
r=154e-6
R = 3.3e-3
compacity = .20
drho = 35.4
p = scontact3.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = -9.81
rhop = 2450
nu = 1.17/1030.
rho = 1030#rhop-drho/compacity
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 1000
#numerical parameters
lcmin = 0.0004 # approx r*100 but should match the mesh size
dt = 5e-2
alpha = 1e-3
epsilon = alpha*lcmin**2 /nu
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
mu = nu*rho
p = scontact3.ParticleProblem()
p.read_vtk(outputdir,2)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 10
outf1 = 100000
#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]
#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Bottom",0,.0),("Bottom",1,0.),("Bottom",2,0.),("X",0,.0),("X",1,0.),("X",2,0.),("Z",0,.0),("Z",1,0.),("Z",2,0.),("Top",3,0.)]
strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
print("Import start!!!")
fluid.import_vtk(outputdir+"/fluid_00002.vtu")
print("Import done!!!")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
tEnd1 = 100
ii = 21
t = 21*5e-5
tic = time.clock()
forces = g*p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
if (ii%5==0 and ii != 0):
fluid.adapt_mesh(5e-3,8e-4,600000)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
for i in range(nsub) :
p.iterate(dt/nsub, forces)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
L = 0.05;
H = 0.6;
P = 0.02;
y = 0;
lc = 0.01;
Point(1) = {-L, H, -P, lc};
Point(2) = {-L, -H, -P, lc};
Point(3) = {L, -H, -P, lc};
Point(4) = {L, H, -P, lc};
Point(5) = {-L, H, P, lc};
Point(6) = {-L, -H, P, lc};
Point(7) = {L, -H, P, lc};
Point(8) = {L, H, P, lc};
Point(9) = {0,H,0,lc};
Point(10) = {0,-H,0,lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line(9) = {5, 1};
Line(10) = {4, 8};
Line(11) = {2, 6};
Line(12) = {3, 7};
Line Loop(1) = {1:4};
Line Loop(2) = {5:8};
Line Loop(3) = {9,-4,10,8};
Line Loop(4) = {2,12,-6,-11};
Line Loop(5) = {12,7,-10,-3};
Line Loop(6) = {1,11,-5,9};
Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Plane Surface(4) = {4};
Plane Surface(5) = {5};
Plane Surface(6) = {6};
Surface Loop(1) = {1,-3,-2,4,1,-5,-2,6};
Volume(1) = {1};
Physical Surface("Z") = {1,2};
Physical Surface("X") = {5,6};
Physical Surface("Top") = {3};
Physical Surface("Bottom") = {4};
Physical Volume("Domain") = {1};
Physical Point("PtFix") = {1};
Point(15) = {0,.52,0};
Field[1] = Attractor;
Field[1].NodesList = {15};
Field[2] = Threshold;
Field[2].DistMax = 0.02;
Field[2].DistMin = 0.01;
Field[2].LcMax = 0.01;
Field[2].LcMin = 0.0008;
Field[2].IField = 1;
Background Field = 2;