Skip to content
Snippets Groups Projects
Commit 7c3cf4b7 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

This commit was manufactured by cvs2svn to create tag 'gmsh_1_63'.

parent 4cacb132
No related branches found
No related tags found
No related merge requests found
// $Id: 3D_Bricks.cpp,v 1.18 2006-01-28 18:44:19 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA.
//
// Please report all bugs and problems to <gmsh@geuz.org>.
#include "Gmsh.h"
#include "Numeric.h"
#include "Mesh.h"
void AddSimplexInGrid(Mesh * m, Simplex * s, int boule_boite)
{
double XminBox = 0., XmaxBox = 0., YminBox = 0., YmaxBox = 0., ZmaxBox = 0., ZminBox = 0.;
int Ix1, Ix2, Iy1, Iy2, Iz1, Iz2;
int i, j, k, index;
Brick Br, *pBrick;
if(!m->Grid.init) {
m->Grid.Bricks =
List_Create(m->Grid.Nx * m->Grid.Ny * m->Grid.Nz, 10, sizeof(Brick));
for(i = 0; i < m->Grid.Nx * m->Grid.Ny * m->Grid.Nz; i++) {
Br.pT = List_Create(2, 2, sizeof(Simplex *));
Br.N = i + 1;
List_Add(m->Grid.Bricks, &Br);
}
m->Grid.init = 1;
}
if(boule_boite == BOITE) {
XminBox = XmaxBox = s->V[0]->Pos.X;
YminBox = YmaxBox = s->V[0]->Pos.Y;
ZminBox = ZmaxBox = s->V[0]->Pos.Z;
for(i = 1; i < 4; i++) {
XminBox = DMIN(XminBox, s->V[i]->Pos.X);
XmaxBox = DMAX(XmaxBox, s->V[i]->Pos.X);
YminBox = DMIN(YminBox, s->V[i]->Pos.Y);
YmaxBox = DMAX(YmaxBox, s->V[i]->Pos.Y);
ZminBox = DMIN(ZminBox, s->V[i]->Pos.Z);
ZmaxBox = DMAX(ZmaxBox, s->V[i]->Pos.Z);
}
}
else if(boule_boite == BOULE) {
XminBox = s->Center.X - s->Radius;
XmaxBox = s->Center.X + s->Radius;
YminBox = s->Center.Y - s->Radius;
YmaxBox = s->Center.Y + s->Radius;
ZminBox = s->Center.Z - s->Radius;
ZmaxBox = s->Center.Z + s->Radius;
}
Ix1 = (int)((double)m->Grid.Nx * (XminBox - m->Grid.min.X) /
(m->Grid.max.X - m->Grid.min.X));
Ix2 = (int)((double)m->Grid.Nx * (XmaxBox - m->Grid.min.X) /
(m->Grid.max.X - m->Grid.min.X));
Iy1 = (int)((double)m->Grid.Ny * (YminBox - m->Grid.min.Y) /
(m->Grid.max.Y - m->Grid.min.Y));
Iy2 = (int)((double)m->Grid.Ny * (YmaxBox - m->Grid.min.Y) /
(m->Grid.max.Y - m->Grid.min.Y));
Iz1 = (int)((double)m->Grid.Nz * (ZminBox - m->Grid.min.Z) /
(m->Grid.max.Z - m->Grid.min.Z));
Iz2 = (int)((double)m->Grid.Nz * (ZmaxBox - m->Grid.min.Z) /
(m->Grid.max.Z - m->Grid.min.Z));
Ix1 = IMAX(Ix1, 0);
Ix2 = IMIN(Ix2, m->Grid.Nx - 1);
Iy1 = IMAX(Iy1, 0);
Iy2 = IMIN(Iy2, m->Grid.Ny - 1);
Iz1 = IMAX(Iz1, 0);
Iz2 = IMIN(Iz2, m->Grid.Nz - 1);
for(i = Ix1; i <= Ix2; i++) {
for(j = Iy1; j <= Iy2; j++) {
for(k = Iz1; k <= Iz2; k++) {
index = i + j * m->Grid.Nx + k * m->Grid.Nx * m->Grid.Ny;
pBrick = (Brick *) List_Pointer(m->Grid.Bricks, index);
List_Add(pBrick->pT, &s);
}
}
}
}
This diff is collapsed.
#ifndef _MGEOM_SEARCH_H
#define _MGEOM_SEARCH_H
#include <vector>
#include <algorithm>
#define MIN(x,y) ((x<y)?(x):(y))
#define MAX(x,y) ((x>y)?(x):(y))
#define TOL 1.e-06
namespace AOMD {
template <class T>
class Brick
{
public:
std::vector<T*> Objects;
Brick(){}
T* operator [] (int i)
{
if(i < 0 || i >= Objects.size())throw i;
T *obj = Objects[i];
return obj;
}
int size(){return Objects.size();}
};
template <class T, class GetBox, class InteriorCheck>
class mGeomSearch
{
Brick<T> *theTable;
int Nx,Ny,Nz;
double Xmin,Xmax,Ymin,Ymax,Zmin,Zmax;
int getBrickId (double X, double Y, double Z)
{
int Ix = (int)((double)Nx * (X-Xmin) / (Xmax-Xmin));
int Iy = (int)((double)Ny * (Y-Ymin) / (Ymax-Ymin));
int Iz = (int)((double)Nz * (Z-Zmin) / (Zmax-Zmin));
Ix = MIN(Ix,Nx-1);
Iy = MIN(Iy,Ny-1);
Iz = MIN(Iz,Nz-1);
if(Ix<0)Ix=0;
if(Iy<0)Iy=0;
if(Iz<0)Iz=0;
int index = Ix + Iy * Nx + Iz * Nx * Ny;
return index;
}
Brick<T> * getBrick (int index)
{
if(index <0 || index >= Nx*Ny*Nz)throw index;
Brick<T> *b = &(theTable[index]);
return b;
}
GetBox getBox;
InteriorCheck interiorCheck;
public:
mGeomSearch (double x1,
double x2,
double y1,
double y2,
double z1,
double z2,
GetBox g,
InteriorCheck i,
int nx = 10,
int ny = 10,
int nz = 10) : Nx(nx), Ny(ny), Nz(nz),getBox(g),interiorCheck(i)
{
Xmin = x1-TOL; Xmax = x2+TOL;
Ymin = y1-TOL; Ymax = y2+TOL;
Zmin = z1-TOL; Zmax = z2+TOL;
theTable = new Brick<T> [Nx*Ny*Nz];
}
~mGeomSearch ()
{
delete [] theTable;
};
Brick<T> *getBrick(double X, double Y, double Z)
{
return (getBrick(getBrickId(X,Y,Z)));
}
T * find (double X[3] , double U[3]);
bool AddObject ( T *);
bool RemoveObject ( T *);
};
template <class T,class GetBox, class InteriorCheck>
T * mGeomSearch<T,GetBox,InteriorCheck> :: find ( double X[3] ,
double U[3] )
{
Brick<T> *b = getBrick (X[0],X[1],X[2]);
if(!b) return 0;
for(int i=0;i<b->size();i++)
{
T * obj = (*b)[i];
if(interiorCheck(obj,X,U))return obj;
}
return 0;
}
template <class T,class GetBox, class InteriorCheck>
bool mGeomSearch<T,GetBox, InteriorCheck> :: AddObject ( T * obj )
{
int Ix1,Ix2,Iy1,Iy2,Iz1,Iz2;
int i,j,k,index;
Brick<T> *pBrick;
/*Template Objects must overload getBox function*/
double min[3];
double max[3];
getBox(obj,min,max);
Ix1 = (int)( (double)Nx * (min[0] - Xmin) /( Xmax - Xmin ));
Ix2 = (int)( (double)Nx * (max[0] - Xmin) /( Xmax - Xmin ));
Iy1 = (int)( (double)Ny * (min[1] - Ymin) /( Ymax - Ymin ));
Iy2 = (int)( (double)Ny * (max[1] - Ymin) /( Ymax - Ymin ));
Iz1 = (int)( (double)Nz * (min[2] - Zmin) /( Zmax - Zmin ));
Iz2 = (int)( (double)Nz * (max[2] - Zmin) /( Zmax - Zmin ));
Ix1 = MAX(Ix1,0);
Ix2 = MIN(Ix2,Nx-1);
Iy1 = MAX(Iy1,0);
Iy2 = MIN(Iy2,Ny-1);
Iz1 = MAX(Iz1,0);
Iz2 = MIN(Iz2,Nz-1);
for(i=Ix1;i<=Ix2;i++){
for(j=Iy1;j<=Iy2;j++){
for(k=Iz1;k<=Iz2;k++){
index = i + j * Nx + k * Nx * Ny;
pBrick = getBrick(index);
pBrick->Objects.push_back(obj);
}
}
}
return true;
}
template <class T,class GetBox, class InteriorCheck>
bool mGeomSearch<T,GetBox, InteriorCheck> :: RemoveObject ( T * obj )
{
int Ix1,Ix2,Iy1,Iy2,Iz1,Iz2;
int i,j,k,index;
Brick<T> *pBrick;
/*Template Objects must overload getBox function*/
double min[3];
double max[3];
getBox(obj,min,max);
Ix1 = (int)( (double)Nx * (min[0] - Xmin) /( Xmax - Xmin ));
Ix2 = (int)( (double)Nx * (max[0] - Xmin) /( Xmax - Xmin ));
Iy1 = (int)( (double)Ny * (min[1] - Ymin) /( Ymax - Ymin ));
Iy2 = (int)( (double)Ny * (max[1] - Ymin) /( Ymax - Ymin ));
Iz1 = (int)( (double)Nz * (min[2] - Zmin) /( Zmax - Zmin ));
Iz2 = (int)( (double)Nz * (max[2] - Zmin) /( Zmax - Zmin ));
Ix1 = MAX(Ix1,0);
Ix2 = MIN(Ix2,Nx-1);
Iy1 = MAX(Iy1,0);
Iy2 = MIN(Iy2,Ny-1);
Iz1 = MAX(Iz1,0);
Iz2 = MIN(Iz2,Nz-1);
for(i=Ix1;i<=Ix2;i++){
for(j=Iy1;j<=Iy2;j++){
for(k=Iz1;k<=Iz2;k++){
index = i + j * Nx + k * Nx * Ny;
pBrick = getBrick(index);
pBrick->Objects.erase ( std::remove (pBrick->Objects.begin(),pBrick->Objects.end(),obj) ,
pBrick->Objects.end () );
}
}
}
return true;
}
} // end of namespace
#undef TOL
#endif
lcar1 = .0275;
length = 0.25;
height = 1.0;
depth = 0.25;
Point(newp) = {length/2,height/2,depth,lcar1}; /* Point 1 */
Point(newp) = {length/2,height/2,0,lcar1}; /* Point 2 */
Point(newp) = {-length/2,height/2,depth,lcar1}; /* Point 3 */
Point(newp) = {-length/2,-height/2,depth,lcar1}; /* Point 4 */
Point(newp) = {length/2,-height/2,depth,lcar1}; /* Point 5 */
Point(newp) = {length/2,-height/2,0,lcar1}; /* Point 6 */
Point(newp) = {-length/2,height/2,0,lcar1}; /* Point 7 */
Point(newp) = {-length/2,-height/2,0,lcar1}; /* Point 8 */
Line(1) = {3,1};
Line(2) = {3,7};
Line(3) = {7,2};
Line(4) = {2,1};
Line(5) = {1,5};
Line(6) = {5,4};
Line(7) = {4,8};
Line(8) = {8,6};
Line(9) = {6,5};
Line(10) = {6,2};
Line(11) = {3,4};
Line(12) = {8,7};
Line Loop(13) = {-6,-5,-1,11};
Plane Surface(14) = {13};
Line Loop(15) = {4,5,-9,10};
Plane Surface(16) = {15};
Line Loop(17) = {-3,-12,8,10};
Plane Surface(18) = {17};
Line Loop(19) = {7,12,-2,11};
Plane Surface(20) = {19};
Line Loop(21) = {-4,-3,-2,1};
Plane Surface(22) = {21};
Line Loop(23) = {8,9,6,7};
Plane Surface(24) = {23};
Surface Loop(25) = {14,24,-18,22,16,-20};
Complex Volume(26) = {25};
n = 21;
/*
sizex = 4;
sizey = 6;
sizez = 4;
*/
sizex = n*length;
sizey = n*height;
sizez = n*depth;
sizex = 9;
sizey = 33;
sizez = 9;
Transfinite Line {2,4,9,7} = sizez;
Transfinite Line {11,5,10,12} = sizey;
Transfinite Line {3,1,6,8} = sizex;
Transfinite Surface {14} = {5,4,3,1};
Transfinite Surface {16} = {6,2,1,5};
Transfinite Surface {18} = {6,2,7,8};
Transfinite Surface {20} = {3,7,8,4};
Transfinite Surface {22} = {3,7,2,1};
Transfinite Surface {24} = {4,8,6,5};
Recombine Surface {14,16,18,20,22,24};
Transfinite Volume {26} = {3,7,2,1,4,8,6,5};
Physical Surface(200) = {14,16,18,20,24};
Physical Volume(100) = {26};
/*
Mesh.use_cut_plane = 1;
Mesh.cut_planea = 0;
Mesh.cut_planeb = 0;
Mesh.cut_planec = 1;
Mesh.cut_planed = -0.125;
*/
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment