// contact_init initialises variables need by contact model
func (o *ElemU) contact_init(edat *inp.ElemData) {

	// vertices on faces with contact
	var contactverts []int
	if len(o.Cell.FaceBcs) > 0 {
		lverts := o.Cell.FaceBcs.GetVerts("contact")
		for _, m := range lverts {
			contactverts = append(contactverts, m)
		}
	}
	o.Nq = len(contactverts)

	// contact flag
	o.HasContact = o.Nq > 0
	if !o.HasContact {
		return
	}

	// vertices on contact face; numbering
	o.ContactId2vid = contactverts
	o.Vid2contactId = utl.IntVals(o.Nu, -1)
	o.Qmap = make([]int, o.Nq)
	for μ, m := range o.ContactId2vid {
		o.Vid2contactId[m] = μ
	}

	// flags
	o.Macaulay, o.βrmp, o.κ = GetContactFaceFlags(edat.Extra)

	// allocate coupling matrices
	o.Kuq = la.MatAlloc(o.Nu, o.Nq)
	o.Kqu = la.MatAlloc(o.Nq, o.Nu)
	o.Kqq = la.MatAlloc(o.Nq, o.Nq)
}

func plot_cone(α float64, preservePrev bool) {
	nu, nv := 11, 21
	l := 1.2
	r := math.Tan(α) * l
	S, T := utl.MeshGrid2D(0, l, 0, 2.0*PI, nu, nv)
	X := la.MatAlloc(nv, nu)
	Y := la.MatAlloc(nv, nu)
	Z := la.MatAlloc(nv, nu)
	u := make([]float64, 3)
	v := make([]float64, 3)
	L := rot_matrix()
	for j := 0; j < nu; j++ {
		for i := 0; i < nv; i++ {
			u[0] = S[i][j] * r * math.Cos(T[i][j])
			u[1] = S[i][j] * r * math.Sin(T[i][j])
			u[2] = S[i][j]
			la.MatVecMul(v, 1, L, u)
			X[i][j], Y[i][j], Z[i][j] = v[0], v[1], v[2]
		}
	}
	pp := 0
	if preservePrev {
		pp = 1
	}
	plt.Wireframe(X, Y, Z, io.Sf("color='b', lw=0.5, preservePrev=%d", pp))
}

// StatTable computes the min, ave, max, and dev of values organised in a table
//  Input:
//   x     -- sample
//   std   -- compute standard deviation (σ) instead of average deviation (adev)
//   withZ -- computes z-matrix as well
//  Convention of indices:
//   0=min  1=ave  2=max  3=dev
//  Output:                        min          ave          max          dev
//                                  ↓            ↓            ↓            ↓
//   x00 x01 x02 x03 x04 x05 → y00=min(x0?) y10=ave(x0?) y20=max(x0?) y30=dev(x0?)
//   x10 x11 x12 x13 x14 x15 → y01=min(x1?) y11=ave(x1?) y21=max(x1?) y31=dev(x1?)
//   x20 x21 x22 x23 x24 x25 → y02=min(x2?) y12=ave(x2?) y22=max(x2?) y32=dev(x2?)
//                                  ↓            ↓            ↓            ↓
//                       min → z00=min(y0?) z01=min(y1?) z02=min(y2?) z03=min(y3?)
//                       ave → z10=ave(y0?) z11=ave(y1?) z12=ave(y2?) z13=ave(y3?)
//                       max → z20=max(y0?) z21=max(y1?) z22=max(y2?) z23=max(y3?)
//                       dev → z30=dev(y0?) z31=dev(y1?) z32=dev(y2?) z33=dev(y3?)
//                                  =            =            =            =
//                       min → z00=min(min) z01=min(ave) z02=min(max) z03=min(dev)
//                       ave → z10=ave(min) z11=ave(ave) z12=ave(max) z13=ave(dev)
//                       max → z20=max(min) z21=max(ave) z22=max(max) z23=max(dev)
//                       dev → z30=dev(min) z31=dev(ave) z32=dev(max) z33=dev(dev)
func StatTable(x [][]float64, std, withZ bool) (y, z [][]float64) {

	// dimensions
	m := len(x)
	if m < 1 {
		return
	}
	n := len(x[0])
	if n < 2 {
		return
	}

	// compute y
	y = la.MatAlloc(4, m)
	for i := 0; i < m; i++ {
		y[0][i], y[1][i], y[2][i], y[3][i] = StatBasic(x[i], std)
	}

	// compute z
	if withZ {
		z = la.MatAlloc(4, 4)
		for i := 0; i < 4; i++ {
			z[0][i], z[1][i], z[2][i], z[3][i] = StatBasic(y[i], std)
		}
	}
	return
}

// PlotDeriv plots derivative dR[i][j][k]du[d] (2D only)
// option =  0 : use CalcBasisAndDerivs
//           1 : use NumericalDeriv
func (o *Nurbs) PlotDeriv(l, d int, args string, npts, option int) {
	lbls := []string{"N\\&dN", "numD"}
	switch o.gnd {
	// curve
	case 1:
	// surface
	case 2:
		xx := la.MatAlloc(npts, npts)
		yy := la.MatAlloc(npts, npts)
		zz := la.MatAlloc(npts, npts)
		du0 := (o.b[0].tmax - o.b[0].tmin) / float64(npts-1)
		du1 := (o.b[1].tmax - o.b[1].tmin) / float64(npts-1)
		drdu := make([]float64, 2)
		for m := 0; m < npts; m++ {
			u0 := o.b[0].tmin + float64(m)*du0
			for n := 0; n < npts; n++ {
				u1 := o.b[1].tmin + float64(n)*du1
				u := []float64{u0, u1}
				x := o.Point(u)
				xx[m][n] = x[0]
				yy[m][n] = x[1]
				switch option {
				case 0:
					o.CalcBasisAndDerivs(u)
					o.GetDerivL(drdu, l)
				case 1:
					o.NumericalDeriv(drdu, u, l)
				}
				zz[m][n] = drdu[d]
			}
		}
		plt.Title(io.Sf("%d,%d:%s", l, d, lbls[option]), "size=10")
		plt.Contour(xx, yy, zz, "fsz=8")
	}
}

// PlotRamp plots the ramp function (contour)
func (o *Plotter) DrawRamp(xmi, xma, ymi, yma float64) {
	if o.Rmpf == nil {
		o.set_empty()
		return
	}
	if o.NptsRmp < 2 {
		o.NptsRmp = 101
	}
	if math.Abs(xma-xmi) < 1e-5 {
		xmi, xma = -0.1, 0.1
	}
	if math.Abs(yma-ymi) < 1e-5 {
		ymi, yma = -0.1, 0.1
	}
	xx := la.MatAlloc(o.NptsRmp, o.NptsRmp)
	yy := la.MatAlloc(o.NptsRmp, o.NptsRmp)
	zz := la.MatAlloc(o.NptsRmp, o.NptsRmp)
	dx := (xma - xmi) / float64(o.NptsRmp-1)
	dy := (yma - ymi) / float64(o.NptsRmp-1)
	for i := 0; i < o.NptsRmp; i++ {
		for j := 0; j < o.NptsRmp; j++ {
			xx[i][j] = xmi + float64(i)*dx
			yy[i][j] = ymi + float64(j)*dy
			zz[i][j] = xx[i][j] - o.Rmpf(xx[i][j]+yy[i][j])
		}
	}
	plt.ContourSimple(xx, yy, zz, "colors=['blue'], linewidths=[2], levels=[0]")
}

func Test_imap(tst *testing.T) {

	//utl.Tsilent = false
	chk.PrintTitle("Test imap")

	for name, shape := range factory {
		gndim := shape.Gndim
		if gndim == 1 {
			continue
		}

		io.Pfyel("--------------------------------- %-6s---------------------------------\n", name)

		// check inverse mapping
		tol := 1e-14
		noise := 0.01
		if name == "tri10" {
			tol = 1e-14
		}
		if shape.FaceNvertsMax > 2 {
			noise = 0.0
		}
		nverts := shape.Nverts
		C := la.MatAlloc(gndim, nverts)
		s := []float64{rand.Float64(), rand.Float64(), rand.Float64()} // scale factors
		la.MatCopy(C, 1.0, shape.NatCoords)
		_ = tol
		io.Pf("nverts:%v\n", nverts)
		io.Pf("gndim:%v\n", gndim)
		for i := 0; i < gndim; i++ {
			for j := 0; j < nverts; j++ {
				C[i][j] *= s[i]
				C[i][j] += noise * rand.Float64() // noise
			}
		}

		r := make([]float64, 3)
		x := make([]float64, 3)
		R := la.MatAlloc(gndim, nverts)

		for j := 0; j < nverts; j++ {
			for i := 0; i < gndim; i++ {
				x[i] = C[i][j]
			}
			err := shape.InvMap(r, x, C)
			io.Pf("r:%v\n", r)
			_ = err
			for i := 0; i < gndim; i++ {
				R[i][j] = r[i]
			}
		}

		chk.Matrix(tst, "checking", tol, R, shape.NatCoords)

		io.PfGreen("OK\n")
	}
}

// PlotBasis plots basis function (2D only)
// option =  0 : use CalcBasis
//           1 : use CalcBasisAndDerivs
//           2 : use RecursiveBasis
func (o *Nurbs) PlotBasis(l int, args string, npts, option int) {
	lbls := []string{"CalcBasis function", "CalcBasisAndDerivs function", "RecursiveBasis function"}
	switch o.gnd {
	// curve
	case 1:
		U := make([]float64, npts)
		S := make([]float64, npts)
		du := (o.b[0].tmax - o.b[0].tmin) / float64(npts-1)
		uvec := []float64{0}
		for m := 0; m < npts; m++ {
			U[m] = o.b[0].tmin + float64(m)*du
			uvec[0] = U[m]
			switch option {
			case 0:
				o.CalcBasis(uvec)
				S[m] = o.GetBasisL(l)
			case 1:
				o.CalcBasisAndDerivs(uvec)
				S[m] = o.GetBasisL(l)
			case 2:
				S[m] = o.RecursiveBasis(uvec, l)
			}
		}
		plt.Plot(U, S, args)
		plt.Gll("$u$", io.Sf("$S_%d$", l), "")
	// surface
	case 2:
		xx := la.MatAlloc(npts, npts)
		yy := la.MatAlloc(npts, npts)
		zz := la.MatAlloc(npts, npts)
		du0 := (o.b[0].tmax - o.b[0].tmin) / float64(npts-1)
		du1 := (o.b[1].tmax - o.b[1].tmin) / float64(npts-1)
		for m := 0; m < npts; m++ {
			u0 := o.b[0].tmin + float64(m)*du0
			for n := 0; n < npts; n++ {
				u1 := o.b[1].tmin + float64(n)*du1
				u := []float64{u0, u1}
				x := o.Point(u)
				xx[m][n] = x[0]
				yy[m][n] = x[1]
				switch option {
				case 0:
					o.CalcBasis(u)
					zz[m][n] = o.GetBasisL(l)
				case 1:
					o.CalcBasisAndDerivs(u)
					zz[m][n] = o.GetBasisL(l)
				case 2:
					zz[m][n] = o.RecursiveBasis(u, l)
				}
			}
		}
		plt.Contour(xx, yy, zz, "fsz=7")
	}
	plt.Title(io.Sf("%s @ %d", lbls[option], l), "size=7")
}

// PlotBasis plots basis function (2D only)
// option =  0 : use CalcBasis
//           1 : use CalcBasisAndDerivs
//           2 : use RecursiveBasis
func (o *Nurbs) PlotBasis(l int, args string, npts, option int) {
	lbls := []string{"Nonly", "N\\&dN", "recN"}
	switch o.gnd {
	// curve
	case 1:
		xx := make([]float64, npts)
		yy := make([]float64, npts)
		du0 := (o.b[0].tmax - o.b[0].tmin) / float64(npts-1)
		for m := 0; m < npts; m++ {
			u0 := o.b[0].tmin + float64(m)*du0
			u := []float64{u0}
			x := o.Point(u)
			xx[m] = x[0]
			switch option {
			case 0:
				o.CalcBasis(u)
				yy[m] = o.GetBasisL(l)
			case 1:
				o.CalcBasisAndDerivs(u)
				yy[m] = o.GetBasisL(l)
			case 2:
				yy[m] = o.RecursiveBasis(u, l)
			}
		}
		plt.Plot(xx, yy, "fsz=8")
	// surface
	case 2:
		xx := la.MatAlloc(npts, npts)
		yy := la.MatAlloc(npts, npts)
		zz := la.MatAlloc(npts, npts)
		du0 := (o.b[0].tmax - o.b[0].tmin) / float64(npts-1)
		du1 := (o.b[1].tmax - o.b[1].tmin) / float64(npts-1)
		for m := 0; m < npts; m++ {
			u0 := o.b[0].tmin + float64(m)*du0
			for n := 0; n < npts; n++ {
				u1 := o.b[1].tmin + float64(n)*du1
				u := []float64{u0, u1}
				x := o.Point(u)
				xx[m][n] = x[0]
				yy[m][n] = x[1]
				switch option {
				case 0:
					o.CalcBasis(u)
					zz[m][n] = o.GetBasisL(l)
				case 1:
					o.CalcBasisAndDerivs(u)
					zz[m][n] = o.GetBasisL(l)
				case 2:
					zz[m][n] = o.RecursiveBasis(u, l)
				}
			}
		}
		plt.Contour(xx, yy, zz, "fsz=8")
	}
	plt.Title(io.Sf("%d:%s", l, lbls[option]), "size=10")
}

// Initialises continues initialisation by generating individuals
//  Optional:  obj  XOR  fcn, nf, ng, nh
func (o *Optimiser) Init(gen Generator_t, obj ObjFunc_t, fcn MinProb_t, nf, ng, nh int) {

	// generic or minimisation problem
	if obj != nil {
		o.ObjFunc = obj
	} else {
		if fcn == nil {
			chk.Panic("either ObjFunc or MinProb must be provided")
		}
		o.Nf, o.Ng, o.Nh, o.MinProb = nf, ng, nh, fcn
		o.ObjFunc = func(sol *Solution, cpu int) {
			o.MinProb(o.F[cpu], o.G[cpu], o.H[cpu], sol.Flt, sol.Int, cpu)
			for i, f := range o.F[cpu] {
				sol.Ova[i] = f
			}
			for i, g := range o.G[cpu] {
				sol.Oor[i] = utl.GtePenalty(g, 0.0, 1) // g[i] ≥ 0
			}
			for i, h := range o.H[cpu] {
				h = math.Abs(h)
				sol.Ova[0] += h
				sol.Oor[o.Ng+i] = utl.GtePenalty(o.EpsH, h, 1) // ϵ ≥ |h[i]|
			}
		}
		o.F = la.MatAlloc(o.Ncpu, o.Nf)
		o.G = la.MatAlloc(o.Ncpu, o.Ng)
		o.H = la.MatAlloc(o.Ncpu, o.Nh)
		o.Nova = o.Nf
		o.Noor = o.Ng + o.Nh
	}

	// calc derived parameters
	o.Generator = gen
	o.CalcDerived()

	// allocate solutions
	o.Solutions = NewSolutions(o.Nsol, &o.Parameters)
	o.Groups = make([]*Group, o.Ncpu)
	for cpu := 0; cpu < o.Ncpu; cpu++ {
		o.Groups[cpu] = new(Group)
		o.Groups[cpu].Init(cpu, o.Ncpu, o.Solutions, &o.Parameters)
	}

	// metrics
	o.Metrics = new(Metrics)
	o.Metrics.Init(o.Nsol, &o.Parameters)

	// auxiliary
	o.tmp = NewSolution(0, 0, &o.Parameters)
	o.cpupairs = utl.IntsAlloc(o.Ncpu/2, 2)
	o.iova0 = -1
	o.ova0 = make([]float64, o.Tf)

	// generate trial solutions
	o.generate_solutions(0)
}

// PlotDeriv plots derivative dR[i][j][k]du[d] (2D only)
// option =  0 : use CalcBasisAndDerivs
//           1 : use NumericalDeriv
func (o *Nurbs) PlotDeriv(l, d int, args string, npts, option int) {
	lbls := []string{"CalcBasisAndDerivs function", "NumericalDeriv function"}
	switch o.gnd {
	// curve
	case 1:
		U := make([]float64, npts)
		G := make([]float64, npts)
		du := (o.b[0].tmax - o.b[0].tmin) / float64(npts-1)
		uvec := []float64{0}
		gvec := []float64{0}
		for m := 0; m < npts; m++ {
			U[m] = o.b[0].tmin + float64(m)*du
			uvec[0] = U[m]
			switch option {
			case 0:
				o.CalcBasisAndDerivs(uvec)
				o.GetDerivL(gvec, l)
			case 1:
				o.NumericalDeriv(gvec, uvec, l)
			}
			G[m] = gvec[0]
		}
		plt.Plot(U, G, args)
		plt.Gll("$u$", io.Sf("$G_%d$", l), "")
	// surface
	case 2:
		xx := la.MatAlloc(npts, npts)
		yy := la.MatAlloc(npts, npts)
		zz := la.MatAlloc(npts, npts)
		du0 := (o.b[0].tmax - o.b[0].tmin) / float64(npts-1)
		du1 := (o.b[1].tmax - o.b[1].tmin) / float64(npts-1)
		drdu := make([]float64, 2)
		for m := 0; m < npts; m++ {
			u0 := o.b[0].tmin + float64(m)*du0
			for n := 0; n < npts; n++ {
				u1 := o.b[1].tmin + float64(n)*du1
				u := []float64{u0, u1}
				x := o.Point(u)
				xx[m][n] = x[0]
				yy[m][n] = x[1]
				switch option {
				case 0:
					o.CalcBasisAndDerivs(u)
					o.GetDerivL(drdu, l)
				case 1:
					o.NumericalDeriv(drdu, u, l)
				}
				zz[m][n] = drdu[d]
			}
		}
		plt.Contour(xx, yy, zz, "fsz=7")
	}
	plt.Title(io.Sf("%s @ %d,%d", lbls[option], l, d), "size=7")
}

// register element
func init() {

	// information allocator
	infogetters["phi"] = func(sim *inp.Simulation, cell *inp.Cell, edat *inp.ElemData) *Info {

		// new info
		var info Info

		nverts := cell.Shp.Nverts
		ykeys := []string{"h"}

		info.Dofs = make([][]string, nverts)
		for m := 0; m < nverts; m++ {
			info.Dofs[m] = ykeys
		}

		info.T1vars = ykeys

		// return information
		return &info
	}

	// element allocator
	eallocators["phi"] = func(sim *inp.Simulation, cell *inp.Cell, edat *inp.ElemData, x [][]float64) Elem {

		// basic data
		var o ElemPhi
		o.Cell = cell
		o.X = x
		o.Nu = o.Cell.Shp.Nverts
		o.Ndim = sim.Ndim

		// integration points
		var err error
		o.IpsElem, o.IpsFace, err = o.Cell.Shp.GetIps(edat.Nip, edat.Nipf)
		if err != nil {
			chk.Panic("cannot allocate integration points of solid element with nip=%d and nipf=%d:\n%v", edat.Nip, edat.Nipf, err)
		}

		// local starred variables
		nip := len(o.IpsElem)
		o.ψs = make([]float64, nip)
		o.PhiSign = make([]float64, nip)

		// scratchpad. computed @ each ip
		o.K = la.MatAlloc(o.Nu, o.Nu)
		o.v_0 = la.MatAlloc(nip, o.Ndim)
		o.reinit = false
		// return new element
		return &o
	}
}

func (o PressCylin) CalcStresses(Pvals []float64, nr int) (R []float64, Sr, St [][]float64) {
	R = utl.LinSpace(o.a, o.b, nr)
	np := len(Pvals)
	Sr = la.MatAlloc(np, nr)
	St = la.MatAlloc(np, nr)
	for i, P := range Pvals {
		c := o.Calc_c(P)
		for j := 0; j < nr; j++ {
			Sr[i][j], St[i][j] = o.Stresses(c, R[j])
		}
	}
	return
}

// Ipoints returns the real coordinates of integration points [nip][ndim]
func (o *ElemUP) Ipoints() (coords [][]float64) {
	coords = la.MatAlloc(len(o.U.IpsElem), o.Ndim)
	for idx, ip := range o.U.IpsElem {
		coords[idx] = o.U.Cell.Shp.IpRealCoords(o.U.X, ip)
	}
	return
}

func Test_2dinteg02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("2dinteg02. bidimensional integral")

	// Γ(1/4, 1)
	gamma_1div4_1 := 0.2462555291934987088744974330686081384629028737277219

	x := utl.LinSpace(0, 1, 11)
	y := utl.LinSpace(0, 1, 11)
	m, n := len(x), len(y)
	f := la.MatAlloc(m, n)
	for i := 0; i < m; i++ {
		for j := 0; j < n; j++ {
			f[i][j] = 8.0 * math.Exp(-math.Pow(x[i], 2)-math.Pow(y[j], 4))
		}
	}
	dx, dy := x[1]-x[0], y[1]-y[0]
	Vt := Trapz2D(dx, dy, f)
	Vs := Simps2D(dx, dy, f)
	Vc := math.Sqrt(math.Pi) * math.Erf(1) * (math.Gamma(1.0/4.0) - gamma_1div4_1)
	io.Pforan("Vt = %v\n", Vt)
	io.Pforan("Vs = %v\n", Vs)
	io.Pfgreen("Vc = %v\n", Vc)
	chk.Scalar(tst, "Vt", 0.0114830435645548, Vt, Vc)
	chk.Scalar(tst, "Vs", 1e-4, Vs, Vc)

}

// Ipoints returns the real coordinates of integration points [nip][ndim]
func (o ElemU) Ipoints() (coords [][]float64) {
	coords = la.MatAlloc(len(o.IpsElem), Global.Ndim)
	for idx, ip := range o.IpsElem {
		coords[idx] = o.Shp.IpRealCoords(o.X, ip)
	}
	return
}

func (o *Plotter) Plot_i_alp(x, y []float64, res []*State, sts [][]float64, last bool) {
	nr := len(res)
	nα := len(res[0].Alp)
	if nα == 0 {
		o.set_empty()
		return
	}
	yy := la.MatAlloc(nα, nr)
	for i := 0; i < nr; i++ {
		x[i] = float64(i)
		for j := 0; j < nα; j++ {
			yy[j][i] = res[i].Alp[j]
		}
	}
	for j := 0; j < nα; j++ {
		lbl := io.Sf("$\\alpha_%d$ "+o.Lbl, j)
		plt.Plot(x, yy[j], io.Sf("'r-', ls='-', clip_on=0, color='%s', marker='%s', label=r'%s'", o.Clr, o.Mrk, lbl))
	}
	if last {
		plt.Gll("$i$", "$\\alpha_k$", "leg_out=1, leg_ncol=4, leg_hlen=2")
		if lims, ok := o.Lims["i,alp"]; ok {
			plt.AxisLims(lims)
		}
	}
}

func (o *Rjoint) debug_print_K() {
	sldNn := o.Sld.Cell.Shp.Nverts
	rodNn := o.Rod.Cell.Shp.Nverts
	K := la.MatAlloc(o.Ny, o.Ny)
	start := o.Sld.Nu
	for i := 0; i < o.Ndim; i++ {
		for m := 0; m < sldNn; m++ {
			r := i + m*o.Ndim
			for j := 0; j < o.Ndim; j++ {
				for n := 0; n < sldNn; n++ {
					c := j + n*o.Ndim
					K[r][c] = o.Kss[r][c]
				}
				for n := 0; n < rodNn; n++ {
					c := j + n*o.Ndim
					K[r][start+c] = o.Ksr[r][c]
					K[start+c][r] = o.Krs[c][r]
				}
			}
		}
	}
	for i := 0; i < o.Ndim; i++ {
		for m := 0; m < rodNn; m++ {
			r := i + m*o.Ndim
			for j := 0; j < o.Ndim; j++ {
				for n := 0; n < rodNn; n++ {
					c := j + n*o.Ndim
					K[start+r][start+c] = o.Krr[r][c]
				}
			}
		}
	}
	la.PrintMat("K", K, "%20.10f", false)
}

// Extrapolator computes the extrapolation matrix for this Shape with a combination of integration points 'ips'
//  Note: E[nverts][nip] must be pre-allocated
func (o *Shape) Extrapolator(E [][]float64, ips []*Ipoint) (err error) {
	la.MatFill(E, 0)
	nip := len(ips)
	N := o.GetShapeMatAtIps(ips)
	if nip < o.Nverts {
		ξ := o.GetNodesNatCoordsMat()
		ξh := o.GetIpsNatCoordsMat(ips)
		ξhi := la.MatAlloc(o.Gndim+1, nip)
		Ni := la.MatAlloc(o.Nverts, nip)
		err = la.MatInvG(Ni, N, 1e-10)
		if err != nil {
			return
		}
		err = la.MatInvG(ξhi, ξh, 1e-10)
		if err != nil {
			return
		}
		ξhξhI := la.MatAlloc(nip, nip) // ξh * inv(ξh)
		for k := 0; k < o.Gndim+1; k++ {
			for j := 0; j < nip; j++ {
				for i := 0; i < nip; i++ {
					ξhξhI[i][j] += ξh[i][k] * ξhi[k][j]
				}
				for i := 0; i < o.Nverts; i++ {
					E[i][j] += ξ[i][k] * ξhi[k][j] // ξ * inv(ξh)
				}
			}
		}
		for i := 0; i < o.Nverts; i++ {
			for j := 0; j < nip; j++ {
				for k := 0; k < nip; k++ {
					I_kj := 0.0
					if j == k {
						I_kj = 1.0
					}
					E[i][j] += Ni[i][k] * (I_kj - ξhξhI[k][j])
				}
			}
		}
	} else {
		err = la.MatInvG(E, N, 1e-10)
		if err != nil {
			return
		}
	}
	return
}

// CheckJ check Jacobian matrix
//  Ouptut: cnd -- condition number (with Frobenius norm)
func (o *NlSolver) CheckJ(x []float64, tol float64, chkJnum, silent bool) (cnd float64, err error) {

	// Jacobian matrix
	var Jmat [][]float64
	if o.useDn {
		Jmat = la.MatAlloc(o.neq, o.neq)
		err = o.JfcnDn(Jmat, x)
		if err != nil {
			return 0, chk.Err(_nls_err5, "dense", err.Error())
		}
	} else {
		if o.numJ {
			err = Jacobian(&o.Jtri, o.Ffcn, x, o.fx, o.w, false)
			if err != nil {
				return 0, chk.Err(_nls_err5, "sparse", err.Error())
			}
		} else {
			err = o.JfcnSp(&o.Jtri, x)
			if err != nil {
				return 0, chk.Err(_nls_err5, "sparse(num)", err.Error())
			}
		}
		Jmat = o.Jtri.ToMatrix(nil).ToDense()
	}
	//la.PrintMat("J", Jmat, "%23g", false)

	// condition number
	cnd, err = la.MatCondG(Jmat, "F", 1e-10)
	if err != nil {
		return cnd, chk.Err(_nls_err6, err.Error())
	}
	if math.IsInf(cnd, 0) || math.IsNaN(cnd) {
		return cnd, chk.Err(_nls_err7, cnd)
	}

	// numerical Jacobian
	if !chkJnum {
		return
	}
	var Jtmp la.Triplet
	ws := make([]float64, o.neq)
	err = o.Ffcn(o.fx, x)
	if err != nil {
		return
	}
	Jtmp.Init(o.neq, o.neq, o.neq*o.neq)
	Jacobian(&Jtmp, o.Ffcn, x, o.fx, ws, false)
	Jnum := Jtmp.ToMatrix(nil).ToDense()
	for i := 0; i < o.neq; i++ {
		for j := 0; j < o.neq; j++ {
			chk.PrintAnaNum(io.Sf("J[%d][%d]", i, j), tol, Jmat[i][j], Jnum[i][j], !silent)
		}
	}
	maxdiff := la.MatMaxDiff(Jmat, Jnum)
	if maxdiff > tol {
		err = chk.Err(_nls_err8, maxdiff)
	}
	return
}

func plot_sphere(preservePrev bool) {
	R := 1.0
	U, V := utl.MeshGrid2D(0, PI/2.0, 0, PI/2.0, NU, NV)
	X, Y, Z := la.MatAlloc(NV, NU), la.MatAlloc(NV, NU), la.MatAlloc(NV, NU)
	for j := 0; j < NU; j++ {
		for i := 0; i < NV; i++ {
			X[i][j] = R * math.Cos(U[i][j]) * math.Sin(V[i][j])
			Y[i][j] = R * math.Sin(U[i][j]) * math.Sin(V[i][j])
			Z[i][j] = R * math.Cos(V[i][j])
		}
	}
	pp := 0
	if preservePrev {
		pp = 1
	}
	plt.Wireframe(X, Y, Z, io.Sf("color='k', lw=0.5, preservePrev=%d", pp))
}