func main() {

	// catch errors
	var tst testing.T
	defer func() {
		if mpi.Rank() == 0 {
			if err := recover(); err != nil {
				io.PfRed("ERROR: %v\n", err)
			}
			if tst.Failed() {
				io.PfRed("test failed\n")
			}
		}
		mpi.Stop(false)
	}()
	mpi.Start(false)

	// start global variables and log
	analysis := fem.NewFEM("data/bh16.sim", "", true, true, false, true, true, 0)

	// run simulation
	err := analysis.Run()
	if err != nil {
		tst.Error("Run failed\n")
		return
	}

	// check
	skipK := true
	tolK := 1e-12
	tolu := 1e-15
	tols := 1e-12
	fem.TestingCompareResultsU(&tst, "data/bh16.sim", "cmp/bh16.cmp", "", tolK, tolu, tols, skipK, true)
}

func main() {

	// catch errors
	var tst testing.T
	defer func() {
		if mpi.Rank() == 0 {
			if err := recover(); err != nil {
				io.PfRed("ERROR: %v\n", err)
			}
			if tst.Failed() {
				io.PfRed("test failed\n")
			}
		}
		mpi.Stop(false)
	}()
	mpi.Start(false)

	// start global variables and log
	analysis := fem.NewFEM("data/p01.sim", "", true, true, false, true, true, 0)

	// run simulation
	err := analysis.Run()
	if err != nil {
		tst.Error("Run failed\n")
		return
	}
}

func check(problem int, tolf, tolg float64) {
	io.Pf("\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> problem %d <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n", problem)
	opt := getfcn(problem)
	x := opt.RptXref
	f := make([]float64, opt.Nf)
	g := make([]float64, opt.Ng)
	h := make([]float64, opt.Nh)
	opt.MinProb(f, g, h, x, nil, 0)
	err := math.Abs(f[0] - opt.RptFref[0])
	io.Pforan("nx=%d nf=%d ng=%d nh=%d\n", opt.Nflt, opt.Nf, opt.Ng, opt.Nh)
	io.Pforan("x = %v\n", x)
	io.Pforan("f = %v  err = %v\n", f, err)
	io.Pforan("g = %v\n", g)
	io.Pforan("h = %v\n", h)
	for i := 0; i < opt.Ng; i++ {
		if g[i] < 0 {
			io.PfRed("unfeasible on g\n")
		}
	}
	for i := 0; i < opt.Nh; i++ {
		if math.Abs(h[i]) > opt.EpsH {
			io.PfRed("unfeasible on h\n")
		}
	}
	if err > tolf {
		io.PfRed("err is too big\n")
	}
}

func main() {

	// catch errors
	var tst testing.T
	defer func() {
		if mpi.Rank() == 0 {
			if err := recover(); err != nil {
				io.PfRed("ERROR: %v\n", err)
			}
			if tst.Failed() {
				io.PfRed("test failed\n")
			}
		}
		mpi.Stop(false)
	}()
	mpi.Start(false)

	// start global variables and log
	if !fem.Start("data/p01.sim", true, true) {
		tst.Error("Start failed\n")
		return
	}

	// make sure to flush log
	defer fem.End()

	// run simulation
	if !fem.Run() {
		tst.Error("Run failed\n")
		return
	}
}

func main() {

	// filename
	filename, fnkey := io.ArgToFilename(0, "sg111", ".sim", true)

	// fem
	if !fem.Start(filename, false, false, false) {
		io.PfRed("Start failed\n")
		return
	}
	dom, sum, ok := fem.AllocSetAndInit(0, true, true)
	if !ok {
		io.PfRed("AllocSetAndInit failed\n")
		return
	}

	// selected node and dof index
	nidx := 1
	didx := 1

	// gofem
	ntout := len(sum.OutTimes)
	uy := make([]float64, ntout)
	for tidx, _ := range sum.OutTimes {

		// read results from file
		if !dom.In(sum, tidx, true) {
			io.PfRed("plot_spo751: cannot read solution\n")
			return
		}

		// collect results for load versus time plot
		nod := dom.Nodes[nidx]
		eq := nod.Dofs[didx].Eq
		uy[tidx] = dom.Sol.Y[eq]

		// check
		if math.Abs(dom.Sol.T-sum.OutTimes[tidx]) > 1e-14 {
			io.PfRed("output times do not match time in solution array\n")
			return
		}
	}

	// plot fem results
	plt.SetForPng(0.8, 400, 200)
	plt.Plot(sum.OutTimes, uy, "'ro-', clip_on=0, label='gofem'")

	// analytical solution
	tAna := utl.LinSpace(0, 5, 101)
	uyAna := make([]float64, len(tAna))
	for i, t := range tAna {
		uyAna[i] = solution_uy(t, 1.0)
	}
	plt.Plot(tAna, uyAna, "'g-', clip_on=0, label='analytical'")

	// save
	plt.Gll("$t$", "$u_y$", "")
	plt.SaveD("/tmp", fnkey+".png")
}

func main() {

	// filename
	filename, fnkey := io.ArgToFilename(0, "rjoint01", ".sim", true)

	// fem
	if !fem.Start(filename, false, false, false) {
		io.PfRed("Start failed\n")
		return
	}
	dom, sum, ok := fem.AllocSetAndInit(0, true, true)
	if !ok {
		io.PfRed("AllocSetAndInit failed\n")
		return
	}

	// rjoint element
	eid := 2
	ele := dom.Elems[eid].(*fem.Rjoint)
	ipd := ele.OutIpsData()

	// load results from file
	n := len(sum.OutTimes)
	mtau0 := make([]float64, n)
	mtau1 := make([]float64, n)
	mtau2 := make([]float64, n)
	ompb0 := make([]float64, n)
	ompb1 := make([]float64, n)
	ompb2 := make([]float64, n)
	for i, _ := range sum.OutTimes {
		if !dom.In(sum, i, true) {
			io.PfRed("cannot read solution\n")
			return
		}
		res0 := ipd[0].Calc(dom.Sol)
		res1 := ipd[1].Calc(dom.Sol)
		res2 := ipd[2].Calc(dom.Sol)
		mtau0[i] = -res0["tau"]
		mtau1[i] = -res1["tau"]
		mtau2[i] = -res2["tau"]
		ompb0[i] = res0["ompb"]
		ompb1[i] = res1["ompb"]
		ompb2[i] = res2["ompb"]
	}

	// plot
	plt.SetForPng(0.8, 400, 200)
	plt.Plot(ompb0, mtau0, "'r-', marker='.', label='p0', clip_on=0")
	plt.Plot(ompb1, mtau1, "'g-', marker='.', label='p1', clip_on=0")
	plt.Plot(ompb2, mtau2, "'b-', marker='.', label='p2', clip_on=0")
	plt.Gll("$\\bar{\\omega}_p$", "$-\\tau$", "")
	plt.SaveD("/tmp", fnkey+".png")
}

// CheckFront0 returns front0 and number of failed/success
func CheckFront0(opt *Optimiser, verbose bool) (nfailed int, front0 []*Solution) {
	front0 = make([]*Solution, 0)
	var nsuccess int
	for _, sol := range opt.Solutions {
		var failed bool
		for _, oor := range sol.Oor {
			if oor > 0 {
				failed = true
				break
			}
		}
		if failed {
			nfailed++
		} else {
			nsuccess++
			if sol.FrontId == 0 {
				front0 = append(front0, sol)
			}
		}
	}
	if verbose {
		if nfailed > 0 {
			io.PfRed("N failed = %d out of %d\n", nfailed, opt.Nsol)
		} else {
			io.PfGreen("N success = %d out of %d\n", nsuccess, opt.Nsol)
		}
		io.PfYel("N front 0 = %d\n", len(front0))
	}
	return
}

// End must be called at the end to flush log file
func End() {
	if err := recover(); err != nil {
		io.PfRed("ERROR: %v\n", err)
	} else {
		fem.End()
	}
}

func main() {

	// dimensions
	lx, ly := 10.0, 3.0

	// define structured mesh data
	var gd gemlab.InData
	gd.Nparts = 4
	gd.Sregs = &gemlab.Sregs{
		Tags: []int{-1},
		Nxs:  []int{10},
		Nys:  []int{3},
		Nzs:  []int{3},
		Points: [][]float64{
			{0, 0}, {lx, 0}, {lx, ly}, {0, ly},
		},
		Conn:  [][]int{{0, 1, 2, 3}},
		Btags: [][]int{{-20, -11, -21, -10}},
	}

	// tag vertices along line (middle vertical line)
	gd.VtagsL = &gemlab.VtagsL{
		Tags: []int{-1},
		Xxa:  [][]float64{{lx / 2.0, 0}},
		Xxb:  [][]float64{{lx / 2.0, ly}},
	}

	fnk := "d2-coarse"
	if err := gemlab.Generate(fnk, &gd); err != nil {
		io.PfRed("%v\n", err.Error())
	}
}

func main() {

	// dimensions
	lx, ly, lz := 10.0, 3.0, 3.0

	// define structured mesh data
	var gd gemlab.InData
	gd.Nparts = 4
	gd.Sregs = &gemlab.Sregs{
		Tags: []int{-1},
		Nxs:  []int{10},
		Nys:  []int{3},
		Nzs:  []int{3},
		Points: [][]float64{
			{0, 0, 0}, {lx, 0, 0}, {lx, ly, 0}, {0, ly, 0},
			{0, 0, lz}, {lx, 0, lz}, {lx, ly, lz}, {0, ly, lz},
		},
		Conn:  [][]int{{0, 1, 2, 3, 4, 5, 6, 7}},
		Btags: [][]int{{-10, -11, -20, -21, -30, -31}},
	}

	fnk := "d3-coarse"
	if err := gemlab.Generate(fnk, &gd); err != nil {
		io.PfRed("%v\n", err.Error())
	}
}

func Test_cxint01(tst *testing.T) {

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

	var ops OpsData
	ops.SetDefault()
	ops.Pc = 1
	ops.Ncuts = 1

	A := []int{1, 2}
	B := []int{-1, -2}
	a := make([]int, len(A))
	b := make([]int, len(A))
	IntCrossover(a, b, A, B, 0, &ops)
	io.Pfred("A = %2d\n", A)
	io.PfRed("B = %2d\n", B)
	io.Pfcyan("a = %2d\n", a)
	io.Pfblue2("b = %2d\n", b)
	chk.Ints(tst, "a", a, []int{1, -2})
	chk.Ints(tst, "b", b, []int{-1, 2})
	io.Pf("\n")

	A = []int{1, 2, 3, 4, 5, 6, 7, 8}
	B = []int{-1, -2, -3, -4, -5, -6, -7, -8}
	a = make([]int, len(A))
	b = make([]int, len(A))
	ops.Cuts = []int{1, 3}
	IntCrossover(a, b, A, B, 0, &ops)
	io.Pfred("A = %2v\n", A)
	io.PfRed("B = %2v\n", B)
	io.Pfcyan("a = %2v\n", a)
	io.Pfblue2("b = %2v\n", b)
	chk.Ints(tst, "a", a, []int{1, -2, -3, 4, 5, 6, 7, 8})
	chk.Ints(tst, "b", b, []int{-1, 2, 3, -4, -5, -6, -7, -8})

	ops.Cuts = []int{5, 7}
	IntCrossover(a, b, A, B, 0, &ops)
	io.Pfred("A = %2v\n", A)
	io.PfRed("B = %2v\n", B)
	io.Pfcyan("a = %2v\n", a)
	io.Pfblue2("b = %2v\n", b)
	chk.Ints(tst, "a", a, []int{1, 2, 3, 4, 5, -6, -7, 8})
	chk.Ints(tst, "b", b, []int{-1, -2, -3, -4, -5, 6, 7, -8})
}

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_out02(tst *testing.T) {

	// finalise analysis process and catch errors
	defer func() {
		if err := recover(); err != nil {
			tst.Fail()
			io.PfRed("ERROR: %v\n", err)
		}
	}()

	// test title
	//verbose()
	chk.PrintTitle("out02")

	// start simulation
	processing := fem.NewFEM("data/twoqua4.sim", "", true, true, false, false, chk.Verbose, 0)

	// run simulation
	err := processing.Run()
	if err != nil {
		tst.Errorf("Run failed:\n%v", err)
		return
	}

	// start post-processing
	Start("data/twoqua4.sim", 0, 0)

	// get second ip coordinates
	xip := Ipoints[1].X
	io.Pfcyan("xip = %v\n", xip)

	// define points
	Define("A", N{-1})
	Define("ips", Along{{xip[0], 0}, {xip[0], 1}})

	// load results
	LoadResults(nil)

	// solution
	var sol ana.CteStressPstrain
	sol.Init(fun.Prms{
		&fun.Prm{N: "qnH", V: -50},
		&fun.Prm{N: "qnV", V: -100},
	})

	// check displacements
	tolu := 1e-15
	x := GetCoords("A")
	ux := GetRes("ux", "A", 0)
	uy := GetRes("uy", "A", 0)
	io.Pforan("ux=%v uy=%v\n", ux, uy)
	for j, t := range Times {
		io.Pfyel("t=%g\n", t)
		sol.CheckDispl(tst, t, []float64{ux[j], uy[j]}, x, tolu)
	}
}

func main() {

	// catch errors
	defer func() {
		if err := recover(); err != nil {
			io.PfRed("ERROR: %v\n", err)
		}
	}()

	// input data
	var mshfn string
	mshfn, fnkey = io.ArgToFilename(0, "data/d2-coarse", ".msh", true)
	io.Pf("\n%s\n", io.ArgsTable(
		"mesh filename", "mshfn", mshfn,
	))

	// read mesh
	msh, err := inp.ReadMsh("", mshfn, 0)
	if err != nil {
		io.PfRed("cannot read mesh:\n%v", err)
		return
	}
	ndim = msh.Ndim
	verts = msh.Verts
	cells = msh.Cells
	dirout = "/tmp/gofem"

	// buffers
	geo := new(bytes.Buffer)
	vtu := new(bytes.Buffer)

	// generate topology
	topology(geo)

	// points data
	pdata_write(vtu)

	// cells data
	cdata_write(vtu)

	// write vtu file
	vtu_write(geo, vtu)
}

func main() {

	// filename
	filename, fnkey := io.ArgToFilename(0, "sg1121", ".sim", true)

	// results
	out.Start(filename, 0, 0)
	out.Define("A", out.N{30})
	out.LoadResults(nil)

	// plot FEM results
	out.Plot("t", "uy", "A", plt.Fmt{C: "k", Ls: "-", L: "gofem"}, -1)

	// old results
	b, err := io.ReadFile("cmp/sg1121gofemold.json")
	if err != nil {
		io.PfRed("cannot read comparison file\n")
		return
	}
	var gofemold struct {
		Time, Uy30 []float64
	}
	err = json.Unmarshal(b, &gofemold)
	if err != nil {
		io.PfRed("cannot unmarshal comparison file\n")
		return
	}

	// mechsys results
	_, res, err := io.ReadTable("cmp/sg1121mechsysN30.cmp")
	if err != nil {
		io.PfRed("cannot read mechsys comparison file\n")
		return
	}

	// save
	plt.SetForPng(0.8, 400, 200)
	out.Draw("/tmp", fnkey+".png", false, func(i, j, n int) {
		plt.Plot(gofemold.Time, gofemold.Uy30, "'r-', lw=2, label='gofemOld'")
		plt.Plot(res["Time"], res["uy"], "'b-', label='mechsys'")
	})
}

func main() {

	// catch errors
	var tst testing.T
	defer func() {
		if mpi.Rank() == 0 {
			if err := recover(); err != nil {
				io.PfRed("ERROR: %v\n", err)
			}
			if tst.Failed() {
				io.PfRed("test failed\n")
			}
		}
		mpi.Stop(false)
	}()
	mpi.Start(false)

	// start global variables and log
	if !fem.Start("data/spo751.sim", true, true) {
		tst.Error("Start failed\n")
		return
	}

	// make sure to flush log
	defer fem.End()

	// run simulation
	if !fem.Run() {
		tst.Error("Run failed\n")
		return
	}

	// check
	skipK := true
	tolK := 1e-17
	tolu := 1e-12
	tols := 1e-14
	fem.TestingCompareResultsU(&tst, "data/spo751.sim", "cmp/spo751.cmp", tolK, tolu, tols, skipK, true)
}

func main() {

	// catch errors
	defer func() {
		if err := recover(); err != nil {
			io.PfRed("ERROR: %v\n", err)
		}
	}()

	// input data
	simfile, _ := io.ArgToFilename(0, "simfile.sim", true)
	zmin := io.ArgToFloat(1, 0.0)
	zmax := io.ArgToFloat(2, 3.0)
	npts := io.ArgToInt(3, 11)
	io.Pf("\n%s\n", io.ArgsTable(
		"simulation filename", "simfile", simfile,
		"min elevation", "zmin", zmin,
		"max elevation", "zmax", zmax,
		"number of points", "npts", npts,
	))

	// sim file
	sim := inp.ReadSim("", simfile, false)
	if sim == nil {
		io.PfRed("cannot read sim file\n")
		return
	}

	// layer
	var lay fem.GeoLayer
	lay.Zmin = zmin
	lay.Zmax = zmax
	lay.Cl = sim.WaterRho0 / sim.WaterBulk
	//if !lay.ReadPorousParameters(sim.Regions[0],
	// TODO

}

func main() {

	// input data
	simfile := "simfile.sim"
	zmin := 0.0
	zmax := 3.0
	npts := 11

	// parse flags
	flag.Parse()
	if len(flag.Args()) > 0 {
		simfile = flag.Arg(0)
	}
	if len(flag.Args()) > 1 {
		zmin = io.Atof(flag.Arg(1))
	}
	if len(flag.Args()) > 2 {
		zmax = io.Atob(flag.Arg(2))
	}
	if len(flag.Args()) > 3 {
		npts = io.Atoi(flag.Arg(3))
	}

	// print input data
	io.Pf("\nInput data\n")
	io.Pf("==========\n")
	io.Pf("  simfile  = %30s // simulation filename\n", simfile)
	io.Pf("  zmin     = %30s // min elevation\n", zmin)
	io.Pf("  zmax     = %30v // max elevation\n", zmax)
	io.Pf("  npts     = %30v // number of points\n", npts)
	io.Pf("\n")

	// sim file
	sim := inp.ReadSim("", simfile, false)
	if sim == nil {
		io.PfRed("cannot read sim file\n")
		return
	}

	// layer
	var lay fem.GeoLayer
	lay.Zmin = zmin
	lay.Zmax = zmax
	lay.Cl = sim.WaterRho0 / sim.WaterBulk
	//if !lay.ReadPorousParameters(sim.Regions[0],
	// TODO

}

// RunFEA runs FE analysis.
func RunFEA(x []float64, cpu int) (lsf, failed float64) {

	// FemData
	o := FEMDATA[cpu]

	// check for NaNs
	defer func() {
		if math.IsNaN(failed) || math.IsNaN(lsf) {
			io.PfRed("x = %+#v\n", x)
			chk.Panic("NaN: failed=%v lsf=%v\n", failed, lsf)
		}
	}()

	// adjust parameters
	for i, v := range o.Vars {
		v.Prm.Set(x[i])
	}
	for _, prm := range o.Sim.AdjDependent {
		if prm.N == "I22" {
			prm.Set(prm.S * prm.Other.V * prm.Other.V)
		}
	}
	o.Dom.RecomputeKM()

	// run
	err := o.Analysis.SolveOneStage(0, true)
	if err != nil {
		failed = 1
		return
	}

	// displacement based limit-state-function
	if o.AwdU > 0 {
		δ := o.Dom.Sol.Y[o.EqsU[0]]
		for i := 1; i < len(o.EqsU); i++ {
			δ = utl.Max(δ, o.Dom.Sol.Y[o.EqsU[i]])
		}
		lsf = o.AwdU - δ
	}
	return
}

func (o *PrincStrainsUp) dbg_plot(eid, ipid int, time float64) {
	if eid == o.DbgEid && ipid == o.DbgIpId {

		// skip if not selected time
		if o.DbgTime > -1 {
			if math.Abs(time-o.DbgTime) > 1e-7 {
				return
			}
		}

		// check Jacobian
		var cnd float64
		io.PfYel("\nchecking Jacobian: eid=%d ipid=%d\n", eid, ipid)
		cnd, err := o.nls.CheckJ(o.x, o.ChkJacTol, true, o.ChkSilent)
		if err != nil {
			io.PfRed("CheckJ failed:\n%v\n", err)
		}
		io.PfYel("after: cnd(J) = %v\n\n", cnd)
		io.Pfcyan("\nN  = %v\n", o.N)
		io.Pfcyan("Nb = %v\n", o.Nb)
		io.Pfcyan("A  = %v\n", o.A)
		io.Pfcyan("h  = %v\n", o.h)
		io.Pfcyan("Mb = %v\n", o.Mb)
		io.Pfcyan("a  = %v\n", o.a)
		io.Pfcyan("b  = %v\n", o.b)
		io.Pfcyan("c  = %v\n", o.c)
		io.Pfcyan("Lσ = %v\n", o.Lσ)
		io.Pforan("Ne = %v\n", o.Ne)
		io.Pforan("Mbe = %v\n", o.Mbe)
		io.Pforan("Fcoef = %v\n\n", o.Fcoef)

		// plot
		var plr Plotter
		plr.SetFig(false, false, 1.5, 400, "/tmp", io.Sf("fig_stress_eid%d_ipid%d", eid, ipid))
		plr.SetModel(o.Mdl)
		plr.PreCor = o.DbgPco
		plr.Plot([]string{"p,q,ys", "oct,ys"}, o.DbgRes, o.DbgSts, true, true)
	}
}