标签：
#ifndef _GAUSS_HPP_
#define _GAUSS_HPP_
#include "config.h"
#include "gauss.h"
#include "geometry/cumath.h"
#include "layout/layout.h"
#include "utils/floatnum.h"
#include "utils/logger.h"
#include <iostream>
void gaussian_surface::surface_block_index::insert(const surface_block &bgs) {
for (const surface_block &old : _list) {
assert(old.block != bgs.block);
if ( !cross(old.shape.box, bgs.shape.box) ) continue;
_index.insert({std::cref(bgs), std::cref(old)});
_index.insert({std::cref(old), std::cref(bgs)});
}
_list.push_back(bgs);
}
int gaussian_surface::surface_block_index::count(const point_t &pt,
const surface_piece &piece) const {
int cnt = 1;
auto blks = sea::val_view(_index.equal_range(piece.from));
for (const surface_block &g : blks) {
assert(g.block != piece.from.block);
space_t d = manhattan_distance(pt, g.shape);
if ( d > SS_ZERO ) continue;
if ( d < SS_ZERO ) return -1;
if ( piece.drc == sdim::END ) {
auto pics = sea::val_view(g.pieces.equal_range(piece.normal));
for (const surface_piece &pic : pics) {
if ( is_on(pt, pic.plane) ) return -1;
}
} else if ( is_positive(piece.drc) ) {
udim ud = abs(piece.drc);
if ( g.shape.box.min.w(ud) == pt.w(ud) ) return -1;
} else {
udim ud = abs(piece.drc);
if ( g.shape.box.max.w(ud) == pt.w(ud) ) return -1;
}
++cnt;
}
return cnt;
}
gauss_point gaussian_surface::get_random_point(Sea::Random &rand) const {
gauss_point result;
result.gtried = 0;
result.sumpdf = 0.0;
while (true) {
++result.gtried;
size_t id = Sea::select(_areas, _areas + _pieces.size(), rand.value());
const surface_piece &pic = _pieces[id];
static constexpr double GN = 64.0;
double p = (floor(rand.value() * GN) + 0.5) / GN;
double q = (floor(rand.value() * GN) + 0.5) / GN;
result.point = select_point(pic.plane, p, q);
if ( pic.from.overflow &&
!contains(_layout.window().space, result.point) ) {
continue;
}
int rpt = _index.count(result.point, pic);
if ( rpt < 0 || (double)rpt * rand.value() > 1.0 ) {
continue;
}
size_t rid = _layout.layers().locate(result.point.z);
diel_t diel;
if ( dielectric_layers::is_diel_id(rid) ) {
diel = _layout.layers().diel_val(rid);
} else {
size_t f = dielectric_layers::face_to_above_diel(rid);
diel = _layout.layers().diel_val(f-1);
diel = (diel + _layout.layers().diel_val(f)) / 2;
}
double pdf = pdftype(diel, pic.from.dist);
result.sumpdf += pdf;
if ( _maxpdf * rand.value() > pdf ) {
continue;
}
result.drc = pic.drc;
result.dist = pic.from.dist;
result.normal = pic.normal;
result.block = pic.from.block;
result.weight = get_diel(diel) / pdf;
break;
}
if ( result.drc == sdim::END ) {
int x = result.normal.x;
int y = result.normal.y;
if ( rand.value() * (double)(abs(x) + abs(y)) < (double)abs(x) ) {
result.drc = x > 0 ? sdim::PX : sdim::NX;
} else {
result.drc = y > 0 ? sdim::PY : sdim::NY;
}
}
return result;
}
void gaussian_surface::initialize(const wire_net *master) {
for (const wire_block *b : master->blocks()) {
generate_bgs(b);
}
for (surface_piece &p : _pieces) {
if ( p.drc != sdim::END ) continue;
normal_t n = p.normal;
n = normal_t(-n.x, -n.y);
p.from.pieces.insert({n, std::cref(p)});
}
for (surface_block &b : _surface_blocks) {
if ( _layout.window().type == circuit_layout::boundary::NEUMANN ) {
b.overflow = !contains(_layout.window().space, b.shape.box);
}
_index.insert(b);
}
size_t n = _pieces.size();
assert(n != 0);
_areas = new double [n];
for (size_t i = 0; i < n; ++i) {
_area += area_of(_pieces[i].plane).val();
_areas[i] = _area;
}
for (size_t i = 0; i < n; ++i) {
_areas[i] /= _area;
}
assert(Sea::FloatNum::equal(_areas[n-1], 1.0));
_areas[n-1] = 1.0;
space_t min_dist = SS_MAX;
diel_t max_diel = ZERO_DIEL;
for (const surface_block &blk : _surface_blocks) {
min_dist = std::min(min_dist, blk.dist);
}
for (size_t m : _layout.layers().diel_ids()) {
max_diel = std::max(max_diel, _layout.layers().diel_val(m));
}
_maxpdf = pdftype(max_diel, min_dist);
}
class update_dist {
private:
const shape_t &_shape;
bool _is_box;
point2d_t hps1[4];
size_t nhp1;
public:
space_t dist;
update_dist(const shape_t &s, space_t d):
_shape(s), _is_box(is_box(s)), nhp1(get_hpoints(s, hps1)), dist(d) {}
void update(const shape_t &other) {
space_t xd = distance(_shape.box, other.box, udim::X);
if ( !(xd < dist) ) return;
space_t yd = distance(_shape.box, other.box, udim::Y);
if ( !(yd < dist) ) return;
space_t zd = distance(_shape.box, other.box, udim::Z);
if ( !(zd < dist) ) return;
if ( _is_box && is_box(other) ) {
dist = std::max(std::max(xd, yd), zd);
return;
}
space_t hd = dist;
for (size_t i = 0; i < nhp1; ++i) {
hd = std::min(hd, horizontal_manhattan_distance(hps1[i], other));
}
point2d_t hps2[4];
size_t nhp2 = get_hpoints(other, hps2);
for (size_t i = 0; i < nhp2; ++i) {
hd = std::min(hd, horizontal_manhattan_distance(hps2[i], _shape));
}
dist = std::max(hd, zd);
}
static bool is_box(const shape_t &s) {
return s.code == shape_code::CUBOID || s.code == shape_code::CYLIND;
}
};
void gaussian_surface::generate_bgs(const wire_block *mblk) {
const shape_t &s1 = mblk->shape();
space_t dist = size_of(s1.box) * 20;
dist = std::min(dist, _layout.wire_width() * 100);
dist = std::min(dist, get_dist_by_layer(s1.box));
dist = std::min(dist, get_dist_by_window(s1.box));
update_dist updater = update_dist(s1, dist);
for (const wire_block *b : _layout.blocks()) {
if ( b->net() == mblk->net() ) continue;
updater.update(b->shape());
}
dist = updater.dist / 2;
if (dist <= SS_ZERO) {
log_error("Error: the master net is touched by other net!\n").nl();
exit(EXIT_FAILURE);
}
switch ( s1.code ) {
case shape_code::CUBOID:
generate_pieces_for_cuboid(mblk, dist);
break;
case shape_code::CYLIND:
generate_pieces_for_cylind(mblk, dist);
break;
case shape_code::NPRISM:
generate_pieces_for_nprism(mblk, dist);
break;
default:
assert(false);
}
}
void gaussian_surface::generate_pieces_for_cuboid(const wire_block *mblk, space_t dist) {
cuboid_t surface = expand(mblk->shape().box, dist);
_surface_blocks.emplace_back(shape_t(surface), dist, mblk);
for (size_t id = 0; id < SIGNED_DIMENSION; ++id) {
sdim sd = (sdim)id;
udim ud = abs(sd);
rectangle_t rc;
rc.base = surface.min;
if ( is_positive(sd) )
rc.base.w(ud, surface.max.w(ud));
rc.nlen = surface.wlen(next(ud));
rc.plen = surface.wlen(prev(ud));
rc.drc = sd;
_pieces.emplace_back(plane_t(rc), sd, _surface_blocks.back());
}
}
void gaussian_surface::generate_pieces_for_cylind(const wire_block *m, space_t d) {
generate_pieces_for_cuboid(m, d);
}
void gaussian_surface::generate_pieces_for_nprism(const wire_block *mblk, space_t dist) {
const nprism_t &p = mblk->shape().nprism();
std::vector<normal_t> nvs(2*p.npts+1);
for (size_t i = 0; i < p.npts; ++i) {
point2d_t s = p.vertex[i], t = p.next_vertex(i);
vector2d_t e = minus(t, s);
auto n = normal_t(sign(e.y), -sign(e.x));
if ( n.x == 0 || n.y == 0 ) {
nvs[2*i+1] = normal_t(n.y + n.x, n.y - n.x);
nvs[2*i+2] = normal_t(n.x - n.y, n.x + n.y);
} else {
nvs[2*i+1] = nvs[2*i+2] = n;
}
}
nvs[0] = nvs[2*p.npts];
point2d_t vs[12];
size_t nv = 0;
for (size_t i = 0; i < p.npts; ++i) {
normal_t n1 = nvs[2*i], n2 = nvs[2*i+1];
vs[nv++] = move(p.vertex[i], dist * n1.x, dist * n1.y);
if ( eq(n1, n2) ) continue;
if ( eq(plus(n1, n2), normal_t(0, 0)) ) {
normal_t n = normal_t(-n1.y, n1.x);
vs[nv++] = move(p.vertex[i], dist * n.x, dist * n.y);
}
vs[nv++] = move(p.vertex[i], dist * n2.x, dist * n2.y);
}
nprism_t surface = nprism_t(nv, p.minz - dist, p.maxz + dist);
assert( nv <= nprism_t::NPTS );
std::copy_n(vs, nv, surface.vertex);
_surface_blocks.emplace_back(shape_t(surface), dist, mblk);
for (size_t i = 1; i < surface.npts - 1; ++i) {
triangle_t tr = triangle_t(surface.minz, surface.vertex[0],
surface.vertex[i], surface.vertex[i+1], sdim::NZ);
_pieces.emplace_back(plane_t(tr), sdim::NZ, _surface_blocks.back());
tr.drc = sdim::PZ;
tr.zpos = surface.maxz;
_pieces.emplace_back(plane_t(tr), sdim::PZ, _surface_blocks.back());
}
for (size_t i = 0; i < surface.npts; ++i) {
// 添加侧面
point2d_t s = surface.vertex[i], t = surface.next_vertex(i);
if ( s.x == t.x ) {
// 两个点的 x 轴相等, 说明面平行于 yOz 平面, 为长方形面
rectangle_t rc;
rc.base = point_t(s.x, std::min(s.y, t.y), surface.minz);
rc.plen = surface.maxz - surface.minz;
rc.nlen = sea::abs(t.y - s.y);
rc.drc = s.y < t.y ? sdim::PX : sdim::NX;
_pieces.emplace_back(plane_t(rc), rc.drc, _surface_blocks.back());
} else if ( s.y == t.y ) {
rectangle_t rc;
rc.base = point_t(std::min(s.x, t.x), s.y, surface.minz);
rc.nlen = surface.maxz - surface.minz;
rc.plen = sea::abs(t.x - s.x);
rc.drc = s.x < t.x ? sdim::NY : sdim::PY;
_pieces.emplace_back(plane_t(rc), rc.drc, _surface_blocks.back());
} else {
incline_t ic;
ic.base = point_t(s.x, s.y, surface.minz);
ic.length = minus(t, s);
ic.height = surface.maxz - surface.minz;
double x = +ic.length.y.val();
double y = -ic.length.x.val();
double l = 32768.0 / sqrt(x*x + y*y);
int nx = (int)round(x * l);
int ny = (int)round(y * l);
_pieces.emplace_back(plane_t(ic), normal_t(nx, ny), _surface_blocks.back());
}
}
}
space_t gaussian_surface::get_dist_by_layer(const cuboid_t &box) const {
const dielectric_layers &layers = _layout.layers();
if ( layers.is_single() ) {
return SS_MAX;
} else if ( _config.layout.multi_dielectric ==
configure::layout_arg::APPROX ) {
return SS_MAX;
}
size_t delta = (_config.layout.multi_dielectric ==
configure::layout_arg::TWO ||
_config.layout.multi_dielectric == configure::layout_arg::APPROX2)
? 2 : 3;
size_t uid = layers.locate(box.min.z);
if ( dielectric_layers::is_diel_id(uid) ) {
uid = dielectric_layers::diel_to_above_face(uid);
}
space_t ud = (box.min.z - layers.face_pos(uid - delta)) * 2;
size_t vid = layers.locate(box.max.z, uid - 1);
if ( dielectric_layers::is_diel_id(vid) ) {
vid = dielectric_layers::diel_to_below_face(vid);
}
space_t vd = (layers.face_pos(vid + delta) - box.max.z) * 2;
return std::min(ud, vd);
}
space_t gaussian_surface::get_dist_by_window(const cuboid_t &box) const {
if ( _layout.window().type == circuit_layout::boundary::NEUMANN ) {
return SS_MAX;
}
const cuboid_t &win = _layout.window().space;
point_t d = minimal(minus(win.max, box.max), minus(box.min, win.min));
return std::min(std::min(d.x, d.y), d.z);
}
double gaussian_surface::pdftype(diel_t diel, space_t dist) const {
switch ( _config.gauss.type ) {
case configure::gauss_arg::UNIFORM: return 1.0;
case configure::gauss_arg::DEFAULT:
case configure::gauss_arg::DIEL: return get_diel(diel);
case configure::gauss_arg::DIST: return 1.0 / ss2d(dist);
case configure::gauss_arg::BOTH: return get_diel(diel) / ss2d(dist);
}
assert(false);
return 0.0;
}
#endif // _GAUSS_HPP_
标签：
原文地址：http://my.oschina.net/donngchao/blog/522661
