-
Notifications
You must be signed in to change notification settings - Fork 149
/
f_cdist.go
86 lines (76 loc) · 2.11 KB
/
f_cdist.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
// Copyright 2016 The Gosl Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package dbf
import (
"math"
"github.com/cpmech/gosl/chk"
)
// Cdist implements the distance from point to a circle (2D) or a sphere (3D)
// where the circle/sphere is implicitly defined by means of F(x) = 0.
// where
// F(x) = sqrt((x-xc) dot (x-xc)) - r
// with r being the radius and xc the coordinates of the centre.
// Thus F > 0 is outside and F < 0 is inside the circle/sphere
type Cdist struct {
xc []float64 // centre; len(xc) = 2 or 3 => must enter "xc", "yc" and "zc" as parameters
r float64 // radius
}
// set allocators database
func init() {
allocators["cdist"] = func() T { return new(Cdist) }
}
// Init initialises the function
func (o *Cdist) Init(prms Params) {
ndim := 2
for _, p := range prms {
if p.N == "zc" {
ndim = 3
break
}
}
o.xc = make([]float64, ndim)
e := prms.Connect(&o.r, "r", "cdist function")
e += prms.Connect(&o.xc[0], "xc", "cdist function")
e += prms.Connect(&o.xc[1], "yc", "cdist function")
if ndim == 3 {
e += prms.Connect(&o.xc[2], "zc", "cdist function")
}
if e != "" {
chk.Panic("%v\n", e)
}
rtol := 1e-10
if o.r < rtol {
chk.Panic("cdist: radius must be greater than %g", rtol)
}
}
// F returns y = F(t, x)
func (o Cdist) F(t float64, x []float64) float64 {
f := (x[0]-o.xc[0])*(x[0]-o.xc[0]) + (x[1]-o.xc[1])*(x[1]-o.xc[1])
if len(x) == 3 {
f += (x[2] - o.xc[2]) * (x[2] - o.xc[2])
}
return math.Sqrt(f) - o.r
}
// G returns ∂y/∂t_cteX = G(t, x)
func (o Cdist) G(t float64, x []float64) float64 {
return 0
}
// H returns ∂²y/∂t²_cteX = H(t, x)
func (o Cdist) H(t float64, x []float64) float64 {
return 0
}
// Grad returns ∇F = ∂y/∂x = Grad(t, x)
func (o Cdist) Grad(v []float64, t float64, x []float64) {
d := (x[0]-o.xc[0])*(x[0]-o.xc[0]) + (x[1]-o.xc[1])*(x[1]-o.xc[1])
if len(x) == 3 {
d += (x[2] - o.xc[2]) * (x[2] - o.xc[2])
}
d = math.Sqrt(d)
v[0] = (x[0] - o.xc[0]) / d
v[1] = (x[1] - o.xc[1]) / d
if len(x) == 3 {
v[2] = (x[2] - o.xc[1]) / d
}
return
}