PolyhedralFace Class Reference

A polygon, for form factor computation. More...

Collaboration diagram for PolyhedralFace:

Public Member Functions
	PolyhedralFace (const std::vector< kvector_t > &_V=std::vector< kvector_t >(), bool _sym_S2=false)
	Sets internal variables for given vertex chain. More...
double	area () const
void	assert_Ci (const PolyhedralFace &other) const
	Throws if deviation from inversion symmetry is detected. Does not check vertices. More...
complex_t	ff (cvector_t q, bool sym_Ci) const
	Returns the contribution ff(q) of this face to the polyhedral form factor. More...
complex_t	ff_2D (cvector_t qpa) const
	Returns the two-dimensional form factor of this face, for use in a prism. More...
complex_t	ff_2D_direct (cvector_t qpa) const
	Two-dimensional form factor, for use in prism, from sum over edge form factors. More...
complex_t	ff_2D_expanded (cvector_t qpa) const
	Two-dimensional form factor, for use in prism, from power series. More...
complex_t	ff_n (int m, cvector_t q) const
	Returns contribution qn*f_n [of order q^(n+1)] from this face to the polyhedral form factor. More...
complex_t	normalProjectionConj (cvector_t q) const
	Returns conj(q)*normal [BasicVector3D::dot is antilinear in 'this' argument]. More...
double	pyramidalVolume () const
double	radius3d () const

Static Public Member Functions
static double	diameter (const std::vector< kvector_t > &V)
	Static method, returns diameter of circle that contains all vertices. More...

Private Member Functions
void	decompose_q (cvector_t q, complex_t &qperp, cvector_t &qpa) const
	Sets qperp and qpa according to argument q and to this polygon's normal. More...
complex_t	edge_sum_ff (cvector_t q, cvector_t qpa, bool sym_Ci) const
	Returns core contribution to analytic 2d form factor. More...
complex_t	expansion (complex_t fac_even, complex_t fac_odd, cvector_t qpa, double abslevel) const
	Returns sum of n>=1 terms of qpa expansion of 2d form factor. More...
complex_t	ff_n_core (int m, cvector_t qpa, complex_t qperp) const
	Returns core contribution to f_n. More...

Private Attributes
std::vector< PolyhedralEdge >	edges
double	m_area
kvector_t	m_normal
	normal vector of this polygon's plane More...
double	m_radius_2d
	radius of enclosing cylinder More...
double	m_radius_3d
	radius of enclosing sphere More...
double	m_rperp
	distance of this polygon's plane from the origin, along 'm_normal' More...
bool	sym_S2
	if true, then edges obtainable by inversion are not provided More...

Static Private Attributes
static int	n_limit_series = 20
static double	qpa_limit_series = 1e-2
	determines when use power series More...

Detailed Description

A polygon, for form factor computation.

Definition at line 60 of file PolyhedralComponents.h.

Constructor & Destructor Documentation

PolyhedralFace()

PolyhedralFace::PolyhedralFace	(	const std::vector< kvector_t > &	V = std::vector<kvector_t>(),
		bool	_sym_S2 = false
	)

Sets internal variables for given vertex chain.

Parameters

V	oriented vertex list
_sym_S2	true if face has a perpedicular two-fold symmetry axis

Definition at line 126 of file PolyhedralComponents.cpp.

                                                                          : sym_S2(_sym_S2)
{
    size_t NV = V.size();
    if (!NV)
        throw std::logic_error("Face with no edges");
    if (NV < 3)
        throw std::logic_error("Face with less than three edges");
 
    // compute radius in 2d and 3d
    m_radius_2d = diameter(V) / 2;
    m_radius_3d = 0;
    for (const kvector_t& v : V)
        m_radius_3d = std::max(m_radius_3d, v.mag());
 
    // Initialize list of 'edges'.
    // Do not create an edge if two vertices are too close to each other.
    // TODO This is implemented in a somewhat sloppy way: we just skip an edge if it would
    //      be too short. This leaves tiny open edges. In a clean implementation, we
    //      rather should merge adjacent vertices before generating edges.
    for (size_t j = 0; j < NV; ++j) {
        size_t jj = (j + 1) % NV;
        if ((V[j] - V[jj]).mag() < 1e-14 * m_radius_2d)
            continue; // distance too short -> skip this edge
        edges.push_back(PolyhedralEdge(V[j], V[jj]));
    }
    size_t NE = edges.size();
    if (NE < 3)
        throw std::invalid_argument("Face has less than three non-vanishing edges");
 
    // compute n_k, rperp
    m_normal = kvector_t();
    for (size_t j = 0; j < NE; ++j) {
        size_t jj = (j + 1) % NE;
        kvector_t ee = edges[j].E().cross(edges[jj].E());
        if (ee.mag2() == 0) {
            throw std::logic_error("Two adjacent edges are parallel");
        }
        m_normal += ee.unit();
    }
    m_normal /= NE;
    m_rperp = 0;
    for (size_t j = 0; j < NV; ++j)
        m_rperp += V[j].dot(m_normal);
    m_rperp /= NV;
    // assert that the vertices lay in a plane
    for (size_t j = 1; j < NV; ++j)
        if (std::abs(V[j].dot(m_normal) - m_rperp) > 1e-14 * m_radius_3d)
            throw std::logic_error("Face is not planar");
    // compute m_area
    m_area = 0;
    for (size_t j = 0; j < NV; ++j) {
        size_t jj = (j + 1) % NV;
        m_area += m_normal.dot(V[j].cross(V[jj])) / 2;
    }
    // only now deal with inversion symmetry
    if (sym_S2) {
        if (NE & 1)
            throw std::logic_error("Odd #edges violates symmetry S2");
        NE /= 2;
        for (size_t j = 0; j < NE; ++j) {
            if (((edges[j].R() - m_rperp * m_normal) + (edges[j + NE].R() - m_rperp * m_normal))
                    .mag()
                > 1e-12 * m_radius_2d)
                throw std::logic_error("Edge centers violate symmetry S2");
            if ((edges[j].E() + edges[j + NE].E()).mag() > 1e-12 * m_radius_2d)
                throw std::logic_error("Edge vectors violate symmetry S2");
        }
        // keep only half of the egdes
        edges.erase(edges.begin() + NE, edges.end());
    }
}

References RealSpace::cross(), diameter(), BasicVector3D< T >::dot(), RealSpace::dot(), edges, m_area, m_normal, m_radius_2d, m_radius_3d, m_rperp, BasicVector3D< T >::mag2(), sym_S2, and BasicVector3D< T >::unit().

Here is the call graph for this function:

Member Function Documentation

area()

double PolyhedralFace::area ( ) const

inline

Definition at line 67 of file PolyhedralComponents.h.

{ return m_area; }

References m_area.

assert_Ci()

void PolyhedralFace::assert_Ci

(

const PolyhedralFace &

other

)

const

Throws if deviation from inversion symmetry is detected. Does not check vertices.

Definition at line 376 of file PolyhedralComponents.cpp.

{
    if (std::abs(m_rperp - other.m_rperp) > 1e-15 * (m_rperp + other.m_rperp))
        throw std::logic_error("Faces with different distance from origin violate symmetry Ci");
    if (std::abs(m_area - other.m_area) > 1e-15 * (m_area + other.m_area))
        throw std::logic_error("Faces with different areas violate symmetry Ci");
    if ((m_normal + other.m_normal).mag() > 1e-14)
        throw std::logic_error("Faces do not have opposite orientation, violating symmetry Ci");
}

References m_area, m_normal, and m_rperp.

decompose_q()

void PolyhedralFace::decompose_q ( cvector_t  q, complex_t &  qperp, cvector_t &  qpa  ) const

private

Sets qperp and qpa according to argument q and to this polygon's normal.

Definition at line 200 of file PolyhedralComponents.cpp.

{
    qperp = m_normal.dot(q);
    qpa = q - qperp * m_normal;
    // improve numeric accuracy:
    qpa -= m_normal.dot(qpa) * m_normal;
    if (qpa.mag() < eps * std::abs(qperp))
        qpa = cvector_t(0., 0., 0.);
}

References BasicVector3D< T >::dot(), m_normal, and BasicVector3D< T >::mag().

Referenced by ff(), and ff_n().

Here is the call graph for this function:

diameter()

double PolyhedralFace::diameter ( const std::vector< kvector_t > &  V )

static

Static method, returns diameter of circle that contains all vertices.

Definition at line 112 of file PolyhedralComponents.cpp.

{
    double diameterFace = 0;
    for (size_t j = 0; j < V.size(); ++j)
        for (size_t jj = j + 1; jj < V.size(); ++jj)
            diameterFace = std::max(diameterFace, (V[j] - V[jj]).mag());
    return diameterFace;
}

Referenced by PolyhedralFace(), and Polyhedron::Polyhedron().

edge_sum_ff()

complex_t PolyhedralFace::edge_sum_ff ( cvector_t  q, cvector_t  qpa, bool  sym_Ci  ) const

private

Returns core contribution to analytic 2d form factor.

Definition at line 282 of file PolyhedralComponents.cpp.

{
    cvector_t prevec = m_normal.cross(qpa); // complex conjugation will take place in .dot
    complex_t sum = 0;
    complex_t vfacsum = 0;
    for (size_t i = 0; i < edges.size(); ++i) {
        const PolyhedralEdge& e = edges[i];
        complex_t qE = e.qE(qpa);
        complex_t qR = e.qR(qpa);
        complex_t Rfac = sym_S2 ? sin(qR) : (sym_Ci ? cos(e.qR(q)) : exp_I(qR));
        complex_t vfac;
        if (sym_S2 || i < edges.size() - 1) {
            vfac = prevec.dot(e.E());
            vfacsum += vfac;
        } else {
            vfac = -vfacsum; // to improve numeric accuracy: qcE_J = - sum_{j=0}^{J-1} qcE_j
        }
        complex_t term = vfac * Math::sinc(qE) * Rfac;
        sum += term;
        //    std::cout << std::scientific << std::showpos << std::setprecision(16)
        //              << "    sum=" << sum << " term=" << term << " vf=" << vfac << " qE=" << qE
        //              << " qR=" << qR << " sinc=" << Math::sinc(qE) << " Rfac=" << Rfac << "\n";
    }
    return sum;
}

References BasicVector3D< T >::cross(), BasicVector3D< T >::dot(), PolyhedralEdge::E(), edges, exp_I(), m_normal, PolyhedralEdge::qE(), PolyhedralEdge::qR(), Math::sinc(), and sym_S2.

Referenced by ff(), and ff_2D_direct().

Here is the call graph for this function:

expansion()

complex_t PolyhedralFace::expansion ( complex_t  fac_even, complex_t  fac_odd, cvector_t  qpa, double  abslevel  ) const

private

Returns sum of n>=1 terms of qpa expansion of 2d form factor.

Definition at line 253 of file PolyhedralComponents.cpp.

{
#ifdef ALGORITHM_DIAGNOSTIC
    polyhedralDiagnosis.algo += 1;
#endif
    complex_t sum = 0;
    complex_t n_fac = I;
    int count_return_condition = 0;
    for (int n = 1; n < n_limit_series; ++n) {
#ifdef ALGORITHM_DIAGNOSTIC
        polyhedralDiagnosis.order = std::max(polyhedralDiagnosis.order, n);
#endif
        complex_t term = n_fac * (n & 1 ? fac_odd : fac_even) * ff_n_core(n, qpa, 0) / qpa.mag2();
        // std::cout << std::setprecision(16) << "    sum=" << sum << " +term=" << term << "\n";
        sum += term;
        if (std::abs(term) <= eps * std::abs(sum) || std::abs(sum) < eps * abslevel)
            ++count_return_condition;
        else
            count_return_condition = 0;
        if (count_return_condition > 2)
            return sum; // regular exit
        n_fac = mul_I(n_fac);
    }
    throw std::runtime_error("Series f(q_pa) not converged");
}

References ff_n_core(), I, BasicVector3D< T >::mag2(), mul_I(), and n_limit_series.

Referenced by ff(), and ff_2D_expanded().

Here is the call graph for this function:

ff()

complex_t PolyhedralFace::ff	(	cvector_t	q,
		bool	sym_Ci
	)		const

Returns the contribution ff(q) of this face to the polyhedral form factor.

Definition at line 310 of file PolyhedralComponents.cpp.

{
    complex_t qperp;
    cvector_t qpa;
    decompose_q(q, qperp, qpa);
    double qpa_red = m_radius_2d * qpa.mag();
    complex_t qr_perp = qperp * m_rperp;
    complex_t ff0 = (sym_Ci ? 2. * I * sin(qr_perp) : exp_I(qr_perp)) * m_area;
    if (qpa_red == 0) {
        return ff0;
    } else if (qpa_red < qpa_limit_series && !sym_S2) {
        // summation of power series
        complex_t fac_even;
        complex_t fac_odd;
        if (sym_Ci) {
            fac_even = 2. * mul_I(sin(qr_perp));
            fac_odd = 2. * cos(qr_perp);
        } else {
            fac_even = exp_I(qr_perp);
            fac_odd = fac_even;
        }
        return ff0 + expansion(fac_even, fac_odd, qpa, std::abs(ff0));
    } else {
        // direct evaluation of analytic formula
        complex_t prefac;
        if (sym_S2)
            prefac = sym_Ci ? -8. * sin(qr_perp) : 4. * mul_I(exp_I(qr_perp));
        else
            prefac = sym_Ci ? 4. : 2. * exp_I(qr_perp);
        // std::cout << "       qrperp=" << qr_perp << " => prefac=" << prefac << "\n";
        return prefac * edge_sum_ff(q, qpa, sym_Ci) / mul_I(qpa.mag2());
    }
}

References decompose_q(), edge_sum_ff(), exp_I(), expansion(), I, m_area, m_radius_2d, m_rperp, BasicVector3D< T >::mag(), BasicVector3D< T >::mag2(), mul_I(), qpa_limit_series, and sym_S2.

Here is the call graph for this function:

ff_2D()

complex_t PolyhedralFace::ff_2D

(

cvector_t

qpa

)

const

Returns the two-dimensional form factor of this face, for use in a prism.

Definition at line 360 of file PolyhedralComponents.cpp.

{
    if (std::abs(qpa.dot(m_normal)) > eps * qpa.mag())
        throw std::logic_error("ff_2D called with perpendicular q component");
    double qpa_red = m_radius_2d * qpa.mag();
    if (qpa_red == 0) {
        return m_area;
    } else if (qpa_red < qpa_limit_series && !sym_S2) {
        return ff_2D_expanded(qpa);
    } else {
        return ff_2D_direct(qpa);
    }
}

References BasicVector3D< T >::dot(), ff_2D_direct(), ff_2D_expanded(), m_area, m_normal, m_radius_2d, BasicVector3D< T >::mag(), qpa_limit_series, and sym_S2.

Here is the call graph for this function:

ff_2D_direct()

complex_t PolyhedralFace::ff_2D_direct

(

cvector_t

qpa

)

const

Two-dimensional form factor, for use in prism, from sum over edge form factors.

Definition at line 353 of file PolyhedralComponents.cpp.

{
    return (sym_S2 ? 4. : 2. / I) * edge_sum_ff(qpa, qpa, false) / qpa.mag2();
}

References edge_sum_ff(), I, BasicVector3D< T >::mag2(), and sym_S2.

Referenced by ff_2D().

Here is the call graph for this function:

ff_2D_expanded()

complex_t PolyhedralFace::ff_2D_expanded

(

cvector_t

qpa

)

const

Two-dimensional form factor, for use in prism, from power series.

Definition at line 346 of file PolyhedralComponents.cpp.

{
    return m_area + expansion(1., 1., qpa, std::abs(m_area));
}

References expansion(), and m_area.

Referenced by ff_2D().

Here is the call graph for this function:

ff_n()

complex_t PolyhedralFace::ff_n	(	int	m,
		cvector_t	q
	)		const

Returns contribution qn*f_n [of order q^(n+1)] from this face to the polyhedral form factor.

Definition at line 231 of file PolyhedralComponents.cpp.

{
    complex_t qn = q.dot(m_normal); // conj(q)*normal (dot is antilinear in 'this' argument)
    if (std::abs(qn) < eps * q.mag())
        return 0.;
    complex_t qperp;
    cvector_t qpa;
    decompose_q(q, qperp, qpa);
    double qpa_mag2 = qpa.mag2();
    if (qpa_mag2 == 0.) {
        return qn * pow(qperp * m_rperp, n) * m_area * ReciprocalFactorialArray[n];
    } else if (sym_S2) {
        return qn * (ff_n_core(n, qpa, qperp) + ff_n_core(n, -qpa, qperp)) / qpa_mag2;
    } else {
        complex_t tmp = ff_n_core(n, qpa, qperp);
        // std::cout << "DBX ff_n " << n << " " << qn << " " << tmp << " " << qpa_mag2 << "\n";
        return qn * tmp / qpa_mag2;
    }
}

References decompose_q(), BasicVector3D< T >::dot(), ff_n_core(), m_area, m_normal, m_rperp, BasicVector3D< T >::mag(), BasicVector3D< T >::mag2(), and sym_S2.

Here is the call graph for this function:

ff_n_core()

complex_t PolyhedralFace::ff_n_core ( int  m, cvector_t  qpa, complex_t  qperp  ) const

private

Returns core contribution to f_n.

Definition at line 212 of file PolyhedralComponents.cpp.

{
    const cvector_t prevec = 2. * m_normal.cross(qpa); // complex conjugation not here but in .dot
    complex_t ret = 0;
    const complex_t qrperp = qperp * m_rperp;
    for (size_t i = 0; i < edges.size(); ++i) {
        const PolyhedralEdge& e = edges[i];
        const complex_t vfac = prevec.dot(e.E());
        const complex_t tmp = e.contrib(m + 1, qpa, qrperp);
        ret += vfac * tmp;
        //     std::cout << std::scientific << std::showpos << std::setprecision(16)
        //               << "DBX ff_n_core " << m << " " << vfac << " " << tmp
        //               << " term=" << vfac * tmp << " sum=" << ret << "\n";
    }
    return ret;
}

References PolyhedralEdge::contrib(), BasicVector3D< T >::cross(), BasicVector3D< T >::dot(), PolyhedralEdge::E(), edges, m_normal, and m_rperp.

Referenced by expansion(), and ff_n().

Here is the call graph for this function:

normalProjectionConj()

complex_t PolyhedralFace::normalProjectionConj ( cvector_t  q ) const

inline

Returns conj(q)*normal [BasicVector3D::dot is antilinear in 'this' argument].

Definition at line 71 of file PolyhedralComponents.h.

{ return q.dot(m_normal); }

References BasicVector3D< T >::dot(), and m_normal.

Here is the call graph for this function:

pyramidalVolume()

double PolyhedralFace::pyramidalVolume ( ) const

inline

Definition at line 68 of file PolyhedralComponents.h.

{ return m_rperp * m_area / 3; }

References m_area, and m_rperp.

radius3d()

double PolyhedralFace::radius3d ( ) const

inline

Definition at line 69 of file PolyhedralComponents.h.

{ return m_radius_3d; }

References m_radius_3d.

Member Data Documentation

edges

std::vector<PolyhedralEdge> PolyhedralFace::edges

private

Definition at line 84 of file PolyhedralComponents.h.

Referenced by PolyhedralFace(), edge_sum_ff(), and ff_n_core().

m_area

double PolyhedralFace::m_area

private

Definition at line 85 of file PolyhedralComponents.h.

Referenced by PolyhedralFace(), area(), assert_Ci(), ff(), ff_2D(), ff_2D_expanded(), ff_n(), and pyramidalVolume().

m_normal

kvector_t PolyhedralFace::m_normal

private

normal vector of this polygon's plane

Definition at line 86 of file PolyhedralComponents.h.

Referenced by PolyhedralFace(), assert_Ci(), decompose_q(), edge_sum_ff(), ff_2D(), ff_n(), ff_n_core(), and normalProjectionConj().

m_radius_2d

double PolyhedralFace::m_radius_2d

private

radius of enclosing cylinder

Definition at line 88 of file PolyhedralComponents.h.

Referenced by PolyhedralFace(), ff(), and ff_2D().

m_radius_3d

double PolyhedralFace::m_radius_3d

private

radius of enclosing sphere

Definition at line 89 of file PolyhedralComponents.h.

Referenced by PolyhedralFace(), and radius3d().

m_rperp

double PolyhedralFace::m_rperp

private

distance of this polygon's plane from the origin, along 'm_normal'

Definition at line 87 of file PolyhedralComponents.h.

Referenced by PolyhedralFace(), assert_Ci(), ff(), ff_n(), ff_n_core(), and pyramidalVolume().

n_limit_series

int PolyhedralFace::n_limit_series = 20

staticprivate

Definition at line 81 of file PolyhedralComponents.h.

Referenced by expansion().

qpa_limit_series

double PolyhedralFace::qpa_limit_series = 1e-2

staticprivate

determines when use power series

Definition at line 80 of file PolyhedralComponents.h.

Referenced by ff(), and ff_2D().

sym_S2

bool PolyhedralFace::sym_S2

private

if true, then edges obtainable by inversion are not provided

Definition at line 83 of file PolyhedralComponents.h.

Referenced by PolyhedralFace(), edge_sum_ff(), ff(), ff_2D(), ff_2D_direct(), and ff_n().

---

The documentation for this class was generated from the following files:

• /home/www/ba/Sample/HardParticle/PolyhedralComponents.h
• /home/www/ba/Sample/HardParticle/PolyhedralComponents.cpp