func debug_print_up_results(d *Domain) {
	io.Pf("\ntime = %23.10f\n", d.Sol.T)
	for _, v := range d.Msh.Verts {
		n := d.Vid2node[v.Id]
		eqpl := n.GetEq("pl")
		equx := n.GetEq("ux")
		equy := n.GetEq("uy")
		var pl, ux, uy float64
		if eqpl >= 0 {
			pl = d.Sol.Y[eqpl]
		}
		if equx >= 0 {
			ux = d.Sol.Y[equx]
		}
		if equy >= 0 {
			uy = d.Sol.Y[equy]
		}
		if math.Abs(pl) < 1e-13 {
			pl = 0
		}
		if math.Abs(ux) < 1e-13 {
			ux = 0
		}
		if math.Abs(uy) < 1e-13 {
			uy = 0
		}
		io.Pf("%3d : pl=%23.10v ux=%23.10f uy=%23.10f\n", v.Id, pl, ux, uy)
	}
}

// PrintDeep3 prints an array of array of array
func PrintDeep3(name string, A [][][]float64) {
	io.Pf("%s = [\n", name)
	for _, a := range A {
		io.Pf("  %v\n", a)
	}
	io.Pf("]\n")
}

func Test_mmMul01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("mmMul01. MatrixMatrix multiplication")

	a := [][]float64{
		{1.0, 2.0, 3.0},
		{0.5, 0.75, 1.5},
	}
	b := [][]float64{
		{0.1, 0.5, 0.5, 0.75},
		{0.2, 2.0, 2.0, 2.0},
		{0.3, 0.5, 0.5, 0.5},
	}
	c := MatAlloc(2, 4)
	MatMul(c, 1, a, b) // c := 1*a*b
	io.Pf("a = %v\n", a)
	io.Pf("b = %v\n", b)
	io.Pf("c = %v\n", c)

	ccor := [][]float64{
		{1.4, 6.0, 6.0, 6.25},
		{0.65, 2.5, 2.5, 2.625},
	}
	chk.Matrix(tst, "c", 1.0e-15, c, ccor)
}

func Test_mtdeb01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("mtdeb01. Deb's mutation")

	var ops OpsData
	ops.SetDefault()
	ops.Pm = 1.0
	ops.Xrange = [][]float64{{-3, 3}, {-4, 4}}
	ops.EnfRange = true

	rnd.Init(0)

	A := []float64{-1, 1}
	io.Pforan("before: A = %v\n", A)
	FltMutationDeb(A, 10, &ops)
	io.Pforan("after:  A = %v\n", A)

	ha0 := rnd.Histogram{Stations: utl.LinSpace(-3, 3, 11)}

	nsamples := 1000
	aa := make([]float64, len(A))
	a0s := make([]float64, nsamples)
	for _, t := range []int{0, 50, 100} {
		for i := 0; i < nsamples; i++ {
			copy(aa, A)
			FltMutationDeb(aa, t, &ops)
			a0s[i] = aa[0]
		}
		ha0.Count(a0s, true)
		io.Pf("\ntime = %d\n", t)
		io.Pf("%s", rnd.TextHist(ha0.GenLabels("%.1f"), ha0.Counts, 60))
	}
}

func Test_groups01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("groups01")

	Init(0)

	ng := 3            // number of groups
	nints := 12        // number of integers
	size := nints / ng // groups size
	ints := utl.IntRange(nints)
	groups := utl.IntsAlloc(ng, size)
	hists := make([]*IntHistogram, ng)
	for i := 0; i < ng; i++ {
		hists[i] = &IntHistogram{Stations: utl.IntRange(nints + 1)}
	}
	IntGetGroups(groups, ints)
	io.Pfcyan("groups = %v\n", groups)
	for i := 0; i < NSAMPLES; i++ {
		IntGetGroups(groups, ints)
		for j := 0; j < ng; j++ {
			check_repeated(groups[j])
			hists[j].Count(groups[j], false)
		}
	}
	for i := 0; i < ng; i++ {
		io.Pf("\n")
		io.Pf(TextHist(hists[i].GenLabels("%d"), hists[i].Counts, 60))
	}
}

func check_constants(tst *testing.T, E, ν, Kcor, Gcor, lcor float64) {

	K, G := Calc_K_from_Enu(E, ν), Calc_G_from_Enu(E, ν)
	l := Calc_l_from_Enu(E, ν)

	io.Pf("E = %v\n", E)
	io.Pf("ν = %v\n", ν)
	io.Pf("K = %v\n", K)
	io.Pf("G = %v\n", G)
	io.Pf("l = %v\n", l)

	chk.Scalar(tst, "KfromEν", 1e-17, K, Kcor)
	chk.Scalar(tst, "GfromEν", 1e-17, G, Gcor)
	chk.Scalar(tst, "lfromEν", 1e-17, l, lcor)

	EfromKG, νfromKG := Calc_E_from_KG(K, G), Calc_nu_from_KG(K, G)
	chk.Scalar(tst, "EfromKG", 1e-17, EfromKG, E)
	chk.Scalar(tst, "νfromKG", 1e-17, νfromKG, ν)

	EfromlG, νfromlG := Calc_E_from_lG(l, G), Calc_nu_from_lG(l, G)
	chk.Scalar(tst, "EfromlG", 1e-17, EfromlG, E)
	chk.Scalar(tst, "νfromlG", 1e-17, νfromlG, ν)

	EfromKν, GfromKν := Calc_E_from_Knu(K, ν), Calc_G_from_Knu(K, ν)
	chk.Scalar(tst, "EfromKν", 1e-17, EfromKν, E)
	chk.Scalar(tst, "GfromKν", 1e-17, GfromKν, G)
}

func (e *Equations) Print() {
	io.Pf("N1 = %v, N2 = %v, N = %v\n", e.N1, e.N2, e.N)
	io.Pf("RF1 (unknown) =\n %v\n", e.RF1)
	io.Pf("FR1 = \n%v\n", e.FR1)
	io.Pf("RF2 (prescribed) =\n %v\n", e.RF2)
	io.Pf("FR2 = \n%v\n", e.FR2)
}

func Test_mylab03a(tst *testing.T) {

	//verbose()
	chk.PrintTitle("mylab03a. ints: min and max. dbls: min and max")

	A := []int{1, 2, 3, -1, -2, 0, 8, -3}
	mi, ma := IntMinMax(A)
	io.Pf("A      = %v\n", A)
	io.Pf("min(A) = %v\n", mi)
	io.Pf("max(A) = %v\n", ma)
	if mi != -3 {
		chk.Panic("min(A) failed")
	}
	if ma != 8 {
		chk.Panic("max(A) failed")
	}

	if Imin(-1, 2) != -1 {
		chk.Panic("Imin(-1,2) failed")
	}
	if Imax(-1, 2) != 2 {
		chk.Panic("Imax(-1,2) failed")
	}
	if Min(-1, 2) != -1.0 {
		chk.Panic("Min(-1,2) failed")
	}
	if Max(-1, 2) != 2.0 {
		chk.Panic("Max(-1,2) failed")
	}
}

func main() {

	// input filename
	_, fnkey := io.ArgToFilename(0, "frame2d", ".sim", true)

	// simple problems
	var opts []*goga.Optimiser
	if fnkey == "simple" {
		io.Pf("\n\n\n")
		//P := []int{1}
		P := utl.IntRange2(1, 19)
		opts = make([]*goga.Optimiser, len(P))
		for i, problem := range P {
			opts[i] = solve_problem(fnkey, problem)
		}
	} else {
		opts = []*goga.Optimiser{solve_problem(fnkey, 0)}
	}
	if opts[0].PlotSet1 {
		return
	}

	if opts[0].Nsamples > 1 {
		io.Pf("\n")
		rpt := goga.NewTexReport(opts)
		rpt.ShowDEC = false
		rpt.Type = 4
		rpt.TextSize = ""
		rpt.Title = "FORM Reliability: " + fnkey
		rpt.Fnkey = "rel-" + fnkey
		rpt.Generate()
	}
}

func Test_invs04(tst *testing.T) {

	//verbose()
	chk.PrintTitle("invs04")

	a := []float64{-10.0, -20.0, -30.0, 4.0 * SQ2, 5.0 * SQ2, 6.0 * SQ2}
	at := Alloc2()
	Man2Ten(at, a)
	io.Pf("a = %v\n", a)
	chk.Matrix(tst, "Man2Ten", 1e-17, at, [][]float64{
		{-10, 4, 6},
		{4, -20, 5},
		{6, 5, -30},
	})

	b := []float64{-88, -77, -55, -3 * SQ2}
	bt := Alloc2()
	Man2Ten(bt, b)
	io.Pf("b = %v\n", b)
	chk.Matrix(tst, "Man2Ten", 1e-17, bt, [][]float64{
		{-88, -3, 0},
		{-3, -77, 0},
		{0, 0, -55},
	})

	ver := chk.Verbose
	run_invs_tests(tst, a, ver)
	run_invs_tests(tst, b, ver)
}

func Test_deep01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("deep01")

	a := Deep3alloc(3, 2, 4)
	for i := 0; i < 3; i++ {
		for j := 0; j < 2; j++ {
			for k := 0; k < 4; k++ {
				if math.Abs(a[i][j][k]) > 1e-17 {
					tst.Errorf("[1;31ma[i][j][k] failed[0m")
				}
			}
		}
	}
	io.Pf("a = %v\n", a)

	b := Deep4alloc(3, 2, 1, 2)
	for i := 0; i < 3; i++ {
		for j := 0; j < 2; j++ {
			for k := 0; k < 1; k++ {
				for l := 0; l < 2; l++ {
					if math.Abs(b[i][j][k][l]) > 1e-17 {
						tst.Errorf("[1;31mb[i][j][k][l] failed[0m")
					}
				}
			}
		}
	}
	io.Pf("b = %v\n", b)
}

func TestAbs(result, expected, absolute_error float64, test_description string) (status int) {
	switch {
	case math.IsNaN(result) || math.IsNaN(expected):
		status = NaN

	case math.IsInf(result, 0) || math.IsInf(expected, 0):
		status = Inf

	case (expected > 0 && expected < DBL_MIN) || (expected < 0 && expected > -DBL_MIN):
		status = NotEqual

	default:
		if math.Abs(result-expected) > absolute_error {
			status = NotEqual
		} else {
			status = Equal
		}
	}
	if test_description != "" {
		io.Pf(test_description)
		switch status {
		case NaN:
			io.Pf(" [1;31mNaN[0m\n  %v observed\n  %v expected.  diff = %v\n", result, expected, result-expected)
		case Inf:
			io.Pf(" [1;31mInf[0m\n  %v observed\n  %v expected.  diff = %v\n", result, expected, result-expected)
		case Equal:
			io.Pf(" [1;32mOk[0m\n  %v observed\n  %v expected.  diff = %v\n", result, expected, result-expected)
		case NotEqual:
			io.Pf(" [1;31mError[0m\n  %v observed\n  %v expected.  diff = %v\n", result, expected, result-expected)
		}
	}
	return
}

// Deep3Deserialize deserializes an array of array of array in column-compressed format
//  Example:
func Deep3Deserialize(I, P []int, S []float64, debug bool) (A [][][]float64) {
	_, nrows, _, ncols := Deep3GetInfo(I, P, S, false)
	A = make([][][]float64, nrows)
	for i := 0; i < nrows; i++ {
		A[i] = make([][]float64, ncols[i])
	}
	iprev := 0 // previous i
	j := 0     // column index
	for l, i := range I {
		nitems := P[l+1] - P[l]
		if i != iprev { // jumped to new column
			j = 0
		}
		if debug {
			io.Pf("l=%v  i=%v  nitems=%v  j=%v\n", l, i, nitems, j)
		}
		for k, p := 0, P[l]; p < P[l+1]; k, p = k+1, p+1 {
			if debug {
				io.Pf("  k=%v  p=%v  s=%v\n", k, p, S[p])
			}
			if k == 0 {
				A[i][j] = make([]float64, nitems)
			}
			A[i][j][k] = S[p]
		}
		iprev = i
		j += 1
	}
	return
}

// StatF computes statistical information corresponding to objective function idxF
func StatF(o *Optimiser, idxF int, verbose bool) (fmin, fave, fmax, fdev float64, F []float64) {
	nsamples := len(o.BestOvas[idxF])
	if nsamples == 0 {
		if verbose {
			io.Pfred("there are no samples for statistical analysis\n")
		}
		return
	}
	F = make([]float64, nsamples)
	if nsamples == 1 {
		F[0] = o.BestOvas[idxF][0]
		fmin, fave, fmax = F[0], F[0], F[0]
		return
	}
	for i, f := range o.BestOvas[idxF] {
		F[i] = f
	}
	fmin, fave, fmax, fdev = rnd.StatBasic(F, true)
	if verbose {
		str := "\n"
		if len(o.RptFref) == o.Nova {
			str = io.Sf(" (%g)\n", o.RptFref[idxF])
		}
		io.Pf("fmin = %g\n", fmin)
		io.Pf("fave = %g"+str, fave)
		io.Pf("fmax = %g\n", fmax)
		io.Pf("fdev = %g\n", fdev)
		o.fix_formatting_data()
		io.Pf(rnd.BuildTextHist(nice(fmin, o.HistNdig)-o.HistDelFmin, nice(fmax, o.HistNdig)+o.HistDelFmax,
			o.HistNsta, F, o.HistFmt, o.HistLen))
	}
	return
}

func checkinput(tst *testing.T, m *Mesh, nverts, ncells int, X [][]float64, vtags, ctags, parts []int, types []string, V [][]int, etags, ftags [][]int) {
	if len(m.Verts) != nverts {
		tst.Errorf("nverts is incorrect: %d != %d", len(m.Verts), nverts)
		return
	}
	if len(m.Cells) != ncells {
		tst.Errorf("ncells is incorrect: %d != %d", len(m.Cells), ncells)
		return
	}
	io.Pfyel("\nvertices:\n")
	for i, v := range m.Verts {
		io.Pf("%+v\n", v)
		chk.Vector(tst, io.Sf("vertex %2d: X", v.Id), 1e-15, v.X, X[v.Id])
		if v.Tag != vtags[i] {
			tst.Errorf("vtag is incorrect: %d != %d", v.Tag, vtags[i])
			return
		}
	}
	io.Pfyel("\ncells:\n")
	for i, c := range m.Cells {
		io.Pf("%+v\n", c)
		if c.Tag != ctags[i] {
			tst.Errorf("ctag is incorrect: %d != %d", c.Tag, ctags[i])
			return
		}
		if c.Part != parts[i] {
			tst.Errorf("part is incorrect: %d != %d", c.Part, parts[i])
			return
		}
		chk.String(tst, types[i], c.Type)
		chk.Ints(tst, io.Sf("cell %2d : V", c.Id), c.V, V[c.Id])
		chk.Ints(tst, io.Sf("cell %2d : edgetags", c.Id), c.EdgeTags, etags[c.Id])
	}
}

func Test_nurbs02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("nurbs02")

	// NURBS
	b := FactoryNurbs2dPlateHole()
	elems := b.Elements()
	nbasis := b.GetElemNumBasis()
	io.Pforan("nbasis = %v\n", nbasis)
	chk.IntAssert(nbasis, 9) // orders := (2,2) => nbasis = (2+1)*(2+1) = 9

	// check basis and elements
	chk.Ints(tst, "elem[0]", elems[0], []int{2, 3, 2, 3})
	chk.Ints(tst, "elem[1]", elems[1], []int{3, 4, 2, 3})
	chk.Ints(tst, "ibasis0", b.IndBasis(elems[0]), []int{0, 1, 2, 4, 5, 6, 8, 9, 10})
	chk.Ints(tst, "ibasis1", b.IndBasis(elems[1]), []int{1, 2, 3, 5, 6, 7, 9, 10, 11})
	chk.IntAssert(b.GetElemNumBasis(), len(b.IndBasis(elems[0])))

	// check derivatives
	b.CheckDerivs(tst, 11, 1e-5, false)

	// refine NURBS
	c := b.KrefineN(2, false)
	elems = c.Elements()
	chk.IntAssert(c.GetElemNumBasis(), len(c.IndBasis(elems[0])))

	// check refined elements
	io.Pf("\n------------ refined -------------\n")
	chk.Ints(tst, "elem[0]", elems[0], []int{2, 3, 2, 3})
	chk.Ints(tst, "elem[1]", elems[1], []int{3, 4, 2, 3})
	chk.Ints(tst, "elem[2]", elems[2], []int{4, 5, 2, 3})
	chk.Ints(tst, "elem[3]", elems[3], []int{5, 6, 2, 3})
	chk.Ints(tst, "elem[4]", elems[4], []int{2, 3, 3, 4})
	chk.Ints(tst, "elem[5]", elems[5], []int{3, 4, 3, 4})
	chk.Ints(tst, "elem[6]", elems[6], []int{4, 5, 3, 4})
	chk.Ints(tst, "elem[7]", elems[7], []int{5, 6, 3, 4})

	// check refined basis
	chk.Ints(tst, "ibasis0", c.IndBasis(elems[0]), []int{0, 1, 2, 6, 7, 8, 12, 13, 14})
	chk.Ints(tst, "ibasis1", c.IndBasis(elems[1]), []int{1, 2, 3, 7, 8, 9, 13, 14, 15})
	chk.Ints(tst, "ibasis2", c.IndBasis(elems[2]), []int{2, 3, 4, 8, 9, 10, 14, 15, 16})
	chk.Ints(tst, "ibasis3", c.IndBasis(elems[3]), []int{3, 4, 5, 9, 10, 11, 15, 16, 17})
	chk.Ints(tst, "ibasis4", c.IndBasis(elems[4]), []int{6, 7, 8, 12, 13, 14, 18, 19, 20})
	chk.Ints(tst, "ibasis5", c.IndBasis(elems[5]), []int{7, 8, 9, 13, 14, 15, 19, 20, 21})
	chk.Ints(tst, "ibasis6", c.IndBasis(elems[6]), []int{8, 9, 10, 14, 15, 16, 20, 21, 22})
	chk.Ints(tst, "ibasis7", c.IndBasis(elems[7]), []int{9, 10, 11, 15, 16, 17, 21, 22, 23})

	// plot
	if chk.Verbose {
		io.Pf("\n------------ plot -------------\n")
		la := 0 + 0*b.n[0]
		lb := 2 + 1*b.n[0]
		PlotNurbs("/tmp/gosl/gm", "nurbs02a.png", b, 41, true, nil)
		PlotNurbsBasis("/tmp/gosl/gm", "nurbs02b.png", b, la, lb)
		PlotNurbsDerivs("/tmp/gosl/gm", "nurbs02c.png", b, la, lb)
		PlotTwoNurbs("/tmp/gosl/gm", "nurbs02d.png", b, c, 41, true, nil)
	}
}

// Print prints the souble-serial-list
func (o DblSlist) Print(fmt string) {
	for i := 0; i < len(o.Ptrs)-1; i++ {
		for j := o.Ptrs[i]; j < o.Ptrs[i+1]; j++ {
			io.Pf(fmt, o.Vals[j])
		}
		io.Pf("\n")
	}
}

func Test_fun12(tst *testing.T) {

	//verbose()
	chk.PrintTitle("fun12. mul")

	cos, err := New("cos", []*Prm{
		&Prm{N: "a", V: 1},
		&Prm{N: "b/pi", V: 2},
		&Prm{N: "c", V: 1},
	})
	if err != nil {
		tst.Errorf("test failed: %v\n", err)
		return
	}

	lin, err := New("lin", []*Prm{
		&Prm{N: "m", V: 0.5},
		&Prm{N: "ts", V: 0},
	})
	if err != nil {
		tst.Errorf("test failed: %v\n", err)
		return
	}

	mul, err := New("mul", []*Prm{
		&Prm{N: "fa", Fcn: cos},
		&Prm{N: "fb", Fcn: lin},
	})
	if err != nil {
		tst.Errorf("test failed: %v\n", err)
		return
	}

	tmin := 0.0
	tmax := 1.0
	xcte := []float64{0, 0, 0}
	//if true {
	if false {
		withG, withH, save, show := true, true, false, true
		plt.Reset()
		PlotT(cos, "/tmp/gosl", "fun-cos-12.png", tmin, tmax, xcte, 41, "", withG, withH, save, show, nil)
		plt.Reset()
		PlotT(lin, "/tmp/gosl", "fun-lin-12.png", tmin, tmax, xcte, 41, "", withG, withH, save, show, nil)
		plt.Reset()
		PlotT(mul, "/tmp/gosl", "fun-mul-12.png", tmin, tmax, xcte, 41, "", withG, withH, save, show, nil)
	}

	sktol := 1e-10
	dtol := 1e-9
	dtol2 := 1e-8
	ver := chk.Verbose
	tskip := []float64{tmin, tmax}
	CheckDerivT(tst, cos, tmin, tmax, xcte, 11, nil, sktol, dtol, dtol2, ver)
	io.Pf("\n")
	CheckDerivT(tst, lin, tmin, tmax, xcte, 11, tskip, sktol, dtol, dtol2, ver)
	io.Pf("\n")
	CheckDerivT(tst, mul, tmin, tmax, xcte, 11, tskip, sktol, dtol, dtol2, ver)
}

func main() {

	// catch errors
	defer func() {
		if err := recover(); err != nil {
			if mpi.Rank() == 0 {
				chk.Verbose = true
				for i := 8; i > 3; i-- {
					chk.CallerInfo(i)
				}
				io.PfRed("ERROR: %v\n", err)
			}
		}
		mpi.Stop(false)
	}()
	mpi.Start(false)

	// default input parameters

	// read input parameters
	fnamepath, _ := io.ArgToFilename(0, "", ".sim", true)
	verbose := io.ArgToBool(1, true)
	erasePrev := io.ArgToBool(2, true)
	saveSummary := io.ArgToBool(3, true)
	allowParallel := io.ArgToBool(4, true)
	alias := io.ArgToString(5, "")

	// message
	if mpi.Rank() == 0 && verbose {
		io.PfWhite("\nGofem v3 -- Go Finite Element Method\n\n")
		io.Pf("Copyright 2015 Dorival Pedroso and Raul Durand. All rights reserved.\n")
		io.Pf("Use of this source code is governed by a BSD-style\n")
		io.Pf("license that can be found in the LICENSE file.\n\n")

		io.Pf("\n%v\n", io.ArgsTable(
			"filename path", "fnamepath", fnamepath,
			"show messages", "verbose", verbose,
			"erase previous results", "erasePrev", erasePrev,
			"save summary", "saveSummary", saveSummary,
			"allow parallel run", "allowParallel", allowParallel,
			"word to add to results", "alias", alias,
		))
	}

	// profiling?
	defer utl.DoProf(false)()

	// analysis data
	readSummary := false
	analysis := fem.NewFEM(fnamepath, alias, erasePrev, saveSummary, readSummary, allowParallel, verbose, 0)

	// run simulation
	err := analysis.Run()
	if err != nil {
		chk.Panic("Run failed:\n%v", err)
	}
}

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")
	}
}