// ReadMsh reads .msh file
func ReadMsh(fnk string) (nurbss []*Nurbs) {

	// read file
	fn := fnk + ".msh"
	buf, err := io.ReadFile(fn)
	if err != nil {
		chk.Panic("ReadMsh cannot read file = '%s', %v'", fn, err)
	}

	// decode
	var dat Data
	err = json.Unmarshal(buf, &dat)
	if err != nil {
		chk.Panic("ReadMsh cannot unmarshal file = '%s', %v'", fn, err)
	}

	// list of vertices
	verts := make([][]float64, len(dat.Verts))
	for i, _ := range dat.Verts {
		verts[i] = make([]float64, 4)
		for j := 0; j < 4; j++ {
			verts[i][j] = dat.Verts[i].C[j]
		}
	}

	// allocate NURBSs
	nurbss = make([]*Nurbs, len(dat.Nurbss))
	for i, v := range dat.Nurbss {
		nurbss[i] = new(Nurbs)
		nurbss[i].Init(v.Gnd, v.Ords, v.Knots)
		nurbss[i].SetControl(verts, v.Ctrls)
	}
	return
}

func PyFile(filename string) {
	b, err := io.ReadFile(filename)
	if err != nil {
		chk.Panic("PyFile failed:\n%v", err)
	}
	io.Ff(&bb, string(b))
}

// Read reads configuration parameters from JSON file
func (o *Parameters) Read(filenamepath string) {
	o.Default()
	b, err := io.ReadFile(filenamepath)
	if err != nil {
		chk.Panic("cannot read parameters file %q", filenamepath)
	}
	err = json.Unmarshal(b, o)
	if err != nil {
		chk.Panic("cannot unmarshal parameters file %q", filenamepath)
	}
	return
}

// NewData allocates and initialises a new FEM data structure
func NewData(fnkey string, cpu int) *FemData {

	// load data
	var o FemData
	buf, err := io.ReadFile("od-" + fnkey + ".json")
	if err != nil {
		chk.Panic("cannot load data:\n%v", err)
	}
	err = json.Unmarshal(buf, &o)
	if err != nil {
		chk.Panic("cannot unmarshal data:\n%v", err)
	}

	// load FEM data
	o.Analysis = fem.NewMain(fnkey+".sim", io.Sf("cpu%d", cpu), false, false, false, false, false, cpu)
	o.Dom = o.Analysis.Domains[0]
	o.Sim = o.Dom.Sim
	o.Vars = o.Dom.Sim.AdjRandom

	// backup dependent variables
	for _, prm := range o.Sim.AdjDependent {
		prm.S = prm.V // copy V into S
	}

	// set stage
	err = o.Analysis.SetStage(0)
	if err != nil {
		chk.Panic("cannot set stage:\n%v", err)
	}

	// equations to track U
	verts := o.Dom.Msh.Verts
	if o.VtagU < 0 {
		verts = o.Dom.Msh.VertTag2verts[o.VtagU]
		if len(verts) < 1 {
			chk.Panic("cannot find vertices with tag = %d\n", o.VtagU)
		}
	}
	o.EqsU = make([]int, len(verts))
	for i, vert := range verts {
		eq := o.Dom.Vid2node[vert.Id].GetEq(o.Ukey)
		if eq < 0 {
			chk.Panic("cannot find equation corresponding to vertex id=%d and ukey=%q", vert.Id, o.Ukey)
		}
		o.EqsU[i] = eq
	}
	return &o
}

func read_pyfem_rod_data(fn string) []PyfemRodData {

	// read file
	b, err := io.ReadFile(fn)
	if err != nil {
		chk.Panic("%v\n", err.Error())
	}

	// decode
	var dat []PyfemRodData
	err = json.Unmarshal(b, &dat)
	if err != nil {
		chk.Panic("%v\n", err.Error())
	}
	return dat
}

// ReadJson reads json file
func (o *Path) ReadJson(ndim int, fname string) (err error) {

	// read file
	b, err := io.ReadFile(fname)
	if err != nil {
		return chk.Err(_path_err11, fname, err)
	}

	// decode
	err = json.Unmarshal(b, o)
	if err != nil {
		return chk.Err(_path_err12, fname, err)
	}

	// set additional information
	return o.Init(ndim)
}

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

// ReadMat reads all materials data from a .mat JSON file
//  Note: returns nil on errors
func ReadMat(dir, fn string) *MatDb {

	// new mat
	var o MatDb

	// read file
	b, err := io.ReadFile(filepath.Join(dir, fn))
	if err != nil {
		return nil
	}

	// decode
	err = json.Unmarshal(b, &o)
	if err != nil {
		return nil
	}

	// results
	return &o
}

// ReadMat reads all materials data from a .mat JSON file
//  Note: returns nil on errors
func ReadMat(dir, fn string) *MatDb {

	// new mat
	var o MatDb

	// read file
	b, err := io.ReadFile(filepath.Join(dir, fn))
	if LogErr(err, "mat: cannot open materials file "+dir+"/"+fn) {
		return nil
	}

	// decode
	if LogErr(json.Unmarshal(b, &o), "mat: cannot unmarshal materials file "+fn) {
		return nil
	}

	// log
	log.Printf("mat: fn=%s nfunctions=%d nmaterials=%d\n", fn, len(o.Functions), len(o.Materials))
	return &o
}

// Read reads mesh
func Read(fn string) (o *Mesh, err error) {

	// new mesh
	o = new(Mesh)

	// read file
	b, err := io.ReadFile(fn)
	if err != nil {
		return
	}

	// decode
	err = json.Unmarshal(b, &o)
	if err != nil {
		return
	}

	// check
	err = o.Check()
	return
}

// ReadConfParams reads configuration parameters from JSON file
func ReadConfParams(filenamepath string) *ConfParams {

	// new params
	var o ConfParams
	o.SetDefault()

	// read file
	b, err := io.ReadFile(filenamepath)
	if err != nil {
		chk.Panic("cannot read parameters file %q", filenamepath)
	}

	// decode
	err = json.Unmarshal(b, &o)
	if err != nil {
		chk.Panic("cannot unmarshal parameters file %q", filenamepath)
	}

	// results
	o.CalcDerived()
	return &o
}

func Test_step02(tst *testing.T) {

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

	buf, err := io.ReadFile("data/beadpanel.step")
	if err != nil {
		tst.Errorf("cannot read file:\n%v", err)
		return
	}
	dat := string(buf)

	var stp STEP
	err = stp.ParseDATA(dat)
	if err != nil {
		tst.Errorf("Parse filed:\n%v", err)
		return
	}

	for _, c := range stp.BScurves {
		io.Pforan("c = %v\n", c)
	}
}

// MatfileNew2Old converts a new mat file to the old mat file format
//  convertsymbols -- convert back symbols with Greek characters to UTF-8
func MatfileNew2Old(dirout string, fnold, fnnew string, convertsymbols bool) {

	// read file
	b, err := io.ReadFile(fnnew)
	if err != nil {
		io.PfRed("cannot open file: %v", err.Error())
		return
	}

	// decode
	var mats_new MatDb
	err = json.Unmarshal(b, &mats_new)
	if err != nil {
		io.PfRed("cannot unmarshal file: %v", err.Error())
		return
	}

	// oldmatdata implements the old data format
	type oldmatdata struct {
		Name  string
		Desc  string
		Model string
		Prms  []string
		Vals  []float64
		Units []string
		Extra string
	}

	// convert
	var mats_old []oldmatdata
	for _, m := range mats_new.Materials {
		var oldmat oldmatdata
		oldmat.Name = m.Name
		oldmat.Desc = m.Desc
		oldmat.Model = m.Model
		oldmat.Extra = m.Extra
		for _, prm := range m.Prms {
			name := prm.N
			if convertsymbols {
				if n, ok := invertedconversiontable[prm.N]; ok {
					name = n
				}
			}
			oldmat.Prms = append(oldmat.Prms, name)
			oldmat.Vals = append(oldmat.Vals, prm.V)
			if len(prm.U) > 0 {
				oldmat.Units = append(oldmat.Units, prm.U)
			}
		}
		mats_old = append(mats_old, oldmat)
	}

	// encode
	buf, err := json.MarshalIndent(mats_old, "", "  ")
	if err != nil {
		return
	}

	// save file
	if dirout == "" {
		io.WriteFileS(fnold, string(buf))
		return
	}
	io.WriteFileSD(dirout, fnold, string(buf))
}

// MatfileOld2New converts an old mat file to new mat file format
//  convertsymbols -- convert symbols with Greek characters to ANSI
func MatfileOld2New(dirout string, fnnew, fnold string, convertsymbols bool) {

	// oldmatdata implements the old data format
	type oldmatdata struct {
		Name  string
		Desc  string
		Model string
		Prms  []string
		Vals  []float64
		Units []string
		Extra string
	}

	// data holder
	var mats_old []oldmatdata

	// read file
	b, err := io.ReadFile(fnold)
	if err != nil {
		chk.Panic("cannot open file: %v", err.Error())
	}

	// decode
	err = json.Unmarshal(b, &mats_old)
	if err != nil {
		chk.Panic("cannot unmarshal file: %v", err.Error())
	}

	// new data holder
	var mats_new MatDb
	for _, m := range mats_old {
		prms := []*fun.Prm{}
		has_units := len(m.Units) == len(m.Prms)
		for i, prm := range m.Prms {
			name := prm
			if convertsymbols {
				if n, ok := conversiontable[prm]; ok {
					name = n
				}
			}
			if has_units {
				prms = append(prms, &fun.Prm{N: name, V: m.Vals[i], U: m.Units[i]})
			} else {
				prms = append(prms, &fun.Prm{N: name, V: m.Vals[i]})
			}
		}
		mat_new := Material{
			Name:  m.Name,
			Desc:  m.Desc,
			Model: m.Model,
			Extra: m.Extra,
			Prms:  prms,
		}
		mats_new.Materials = append(mats_new.Materials, &mat_new)
	}

	// save file
	if dirout == "" {
		io.WriteFileS(fnnew, mats_new.String())
		return
	}
	io.WriteFileSD(dirout, fnnew, mats_new.String())
}

// NewData allocates and initialises a new FEM data structure
func NewData(filename, fnkey string, cpu int) *FemData {

	// load opt data
	var o FemData
	buf, err := io.ReadFile("od-" + fnkey + ".json")
	if err != nil {
		chk.Panic("cannot load opt data:\n%v", err)
	}
	err = json.Unmarshal(buf, &o.Opt)
	if err != nil {
		chk.Panic("cannot unmarshal opt data:\n%v", err)
	}
	o.Opt.CalcDerived()

	// load FEM data
	o.Analysis = fem.NewFEM(filename, io.Sf("cpu%d", cpu), false, false, false, false, false, cpu)
	o.Reg = o.Analysis.Sim.Regions[0]
	o.Dom = o.Analysis.Domains[0]

	// required vertices
	for _, vtag := range o.Opt.ReqVtags {
		verts := o.Dom.Msh.VertTag2verts[vtag]
		for _, vert := range verts {
			o.ReqVids = append(o.ReqVids, vert.Id)
		}
	}

	// groups
	o.Ncells = len(o.Dom.Msh.Cells)
	o.Cid2xid = make([]int, o.Ncells)
	if o.Opt.Groups {
		o.Nareas = len(o.Dom.Msh.CellTag2cells)
		for tag, cells := range o.Dom.Msh.CellTag2cells {
			xid := -tag - 1
			for _, cell := range cells {
				o.Cid2xid[cell.Id] = xid
			}
		}
	} else { // all cells
		o.Nareas = o.Ncells
		for i, _ := range o.Dom.Msh.Cells {
			o.Cid2xid[i] = i
		}
	}

	// vertex to track deflection
	o.VidU = -1 // allvertices
	if o.Opt.VtagU < 0 {
		verts := o.Dom.Msh.VertTag2verts[o.Opt.VtagU]
		if len(verts) != 1 {
			chk.Panic("cannot set vertex to track deflection: verts = %v", verts)
		}
		o.VidU = verts[0].Id
	}

	// compute max weight
	o.Analysis.SetStage(0)
	for _, elem := range o.Dom.Elems {
		ele := elem.(*fem.ElastRod)
		o.MaxWeight += ele.Mdl.Rho * ele.Mdl.A * ele.L
	}
	return &o
}

// ReadMsh reads a mesh for FE analyses
//  Note: returns nil on errors
func ReadMsh(dir, fn string) *Mesh {

	// new mesh
	var o Mesh

	// read file
	o.FnamePath = filepath.Join(dir, fn)
	b, err := io.ReadFile(o.FnamePath)
	if LogErr(err, "msh: cannot open mesh file "+o.FnamePath) {
		return nil
	}

	// decode
	if LogErr(json.Unmarshal(b, &o), "msh: cannot unmarshal mesh file "+fn+"\n") {
		return nil
	}

	// check
	if LogErrCond(len(o.Verts) < 2, "msh: mesh must have at least 2 vertices and 1 cell") {
		return nil
	}
	if LogErrCond(len(o.Cells) < 1, "msh: mesh must have at least 2 vertices and 1 cell") {
		return nil
	}

	// vertex related derived data
	o.Ndim = 2
	o.Xmin = o.Verts[0].C[0]
	o.Ymin = o.Verts[0].C[1]
	if len(o.Verts[0].C) > 2 {
		o.Zmin = o.Verts[0].C[2]
	}
	o.Xmax = o.Xmin
	o.Ymax = o.Ymin
	o.Zmax = o.Zmin
	o.VertTag2verts = make(map[int][]*Vert)
	for i, v := range o.Verts {

		// check vertex id
		if LogErrCond(v.Id != i, "msh: vertices must be sequentially numbered. %d != %d\n", v.Id, i) {
			return nil
		}

		// ndim
		nd := len(v.C)
		if LogErrCond(nd < 2 || nd > 4, "msh: ndim must be 2 or 3\n") {
			return nil
		}
		if nd == 3 {
			if math.Abs(v.C[2]) > Ztol {
				o.Ndim = 3
			}
		}

		// tags
		if v.Tag < 0 {
			verts := o.VertTag2verts[v.Tag]
			o.VertTag2verts[v.Tag] = append(verts, v)
		}

		// limits
		o.Xmin = min(o.Xmin, v.C[0])
		o.Xmax = max(o.Xmax, v.C[0])
		o.Ymin = min(o.Ymin, v.C[1])
		o.Ymax = max(o.Ymax, v.C[1])
		if nd > 2 {
			o.Zmin = min(o.Zmin, v.C[2])
			o.Zmax = max(o.Zmax, v.C[2])
		}
	}

	// derived data
	o.CellTag2cells = make(map[int][]*Cell)
	o.FaceTag2cells = make(map[int][]CellFaceId)
	o.FaceTag2verts = make(map[int][]int)
	o.SeamTag2cells = make(map[int][]CellSeamId)
	o.Ctype2cells = make(map[string][]*Cell)
	o.Part2cells = make(map[int][]*Cell)
	for i, c := range o.Cells {

		// check id and tag
		if LogErrCond(c.Id != i, "msh: cells must be sequentially numbered. %d != %d\n", c.Id, i) {
			return nil
		}
		if LogErrCond(c.Tag >= 0, "msh: cell tags must be negative\n") {
			return nil
		}

		// face tags
		cells := o.CellTag2cells[c.Tag]
		o.CellTag2cells[c.Tag] = append(cells, c)
		for i, ftag := range c.FTags {
			if ftag < 0 {
				pairs := o.FaceTag2cells[ftag]
				o.FaceTag2cells[ftag] = append(pairs, CellFaceId{c, i})
				for _, l := range shp.GetFaceLocalVerts(c.Type, i) {
					utl.IntIntsMapAppend(&o.FaceTag2verts, ftag, o.Verts[c.Verts[l]].Id)
				}
//.........这里部分代码省略.........

// testing_compare_results_u compares results with u-formulation
func TestingCompareResultsU(tst *testing.T, simfname, cmpfname string, tolK, tolu, tols float64, skipK, verbose bool) {

	// only root can run this test
	if !Global.Root {
		return
	}

	// read summary
	sum := ReadSum(Global.Dirout, Global.Fnkey)
	if sum == nil {
		tst.Error("cannot read summary file for simulation=%q\n", simfname)
		return
	}

	// allocate domain
	distr := false
	d := NewDomain(Global.Sim.Regions[0], distr)
	if !d.SetStage(0, Global.Sim.Stages[0], distr) {
		tst.Errorf("TestingCompareResultsU: SetStage failed\n")
		return
	}

	// read file
	buf, err := io.ReadFile(cmpfname)
	if err != nil {
		tst.Errorf("TestingCompareResultsU: ReadFile failed\n")
		return
	}

	// unmarshal json
	var cmp_set T_results_set
	err = json.Unmarshal(buf, &cmp_set)
	if err != nil {
		tst.Errorf("TestingCompareResultsU: Unmarshal failed\n")
		return
	}

	// run comparisons
	dmult := 1.0
	for idx, cmp := range cmp_set {

		// displacements multiplier
		if idx == 0 && math.Abs(cmp.DispMult) > 1e-10 {
			dmult = cmp.DispMult
		}

		// time index
		tidx := idx + 1
		if verbose {
			io.PfYel("\n\ntidx = %d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n", tidx)
		}

		// load gofem results
		if !d.In(sum, tidx, true) {
			tst.Errorf("TestingCompareResultsU: reading of results failed\n")
			return
		}
		if verbose {
			io.Pfyel("time = %v\n", d.Sol.T)
		}

		// check K matrices
		if !skipK {
			if verbose {
				io.Pfgreen(". . . checking K matrices . . .\n")
			}
			for eid, Ksg := range cmp.Kmats {
				if e, ok := d.Elems[eid].(*ElemU); ok {
					if !e.AddToKb(d.Kb, d.Sol, true) {
						tst.Errorf("TestingCompareResultsU: AddToKb failed\n")
						return
					}
					chk.Matrix(tst, io.Sf("K%d", eid), tolK, e.K, Ksg)
				}
			}
		}

		// check displacements
		if verbose {
			io.Pfgreen(". . . checking displacements . . .\n")
		}
		for nid, usg := range cmp.Disp {
			ix := d.Vid2node[nid].Dofs[0].Eq
			iy := d.Vid2node[nid].Dofs[1].Eq
			chk.AnaNum(tst, "ux", tolu, d.Sol.Y[ix], usg[0]*dmult, verbose)
			chk.AnaNum(tst, "uy", tolu, d.Sol.Y[iy], usg[1]*dmult, verbose)
			if len(usg) == 3 {
				iz := d.Vid2node[nid].Dofs[2].Eq
				chk.AnaNum(tst, "uz", tolu, d.Sol.Y[iz], usg[2]*dmult, verbose)
			}
		}

		// check stresses
		if true {
			if verbose {
				io.Pfgreen(". . . checking stresses . . .\n")
			}
			for eid, sig := range cmp.Sigmas {
				if verbose {
//.........这里部分代码省略.........

// ReadMsh reads a mesh for FE analyses
//  Note: returns nil on errors
func ReadMsh(dir, fn string, goroutineId int) (o *Mesh, err error) {

	// new mesh
	o = new(Mesh)

	// read file
	o.FnamePath = filepath.Join(dir, fn)
	b, err := io.ReadFile(o.FnamePath)
	if err != nil {
		return
	}

	// decode
	err = json.Unmarshal(b, &o)
	if err != nil {
		return
	}

	// check
	if len(o.Verts) < 2 {
		err = chk.Err("at least 2 vertices are required in mesh\n")
		return
	}
	if len(o.Cells) < 1 {
		err = chk.Err("at least 1 cell is required in mesh\n")
		return
	}

	// variables for NURBS
	var controlpts [][]float64
	has_nurbs := false
	if len(o.Nurbss) > 0 {
		controlpts = make([][]float64, len(o.Verts))
		has_nurbs = true
	}

	// vertex related derived data
	o.Ndim = 2
	o.Xmin = o.Verts[0].C[0]
	o.Ymin = o.Verts[0].C[1]
	if len(o.Verts[0].C) > 2 {
		o.Zmin = o.Verts[0].C[2]
	}
	o.Xmax = o.Xmin
	o.Ymax = o.Ymin
	o.Zmax = o.Zmin
	o.VertTag2verts = make(map[int][]*Vert)
	for i, v := range o.Verts {

		// check vertex id
		if v.Id != i {
			err = chk.Err("vertices ids must coincide with order in \"verts\" list. %d != %d\n", v.Id, i)
			return
		}

		// ndim
		nd := len(v.C)
		if nd < 2 || nd > 4 {
			err = chk.Err("number of space dimensions must be 2, 3 or 4 (NURBS). %d is invalid\n", nd)
			return
		}
		if nd == 3 {
			if math.Abs(v.C[2]) > Ztol {
				o.Ndim = 3
			}
		}

		// tags
		if v.Tag < 0 {
			verts := o.VertTag2verts[v.Tag]
			o.VertTag2verts[v.Tag] = append(verts, v)
		}

		// limits
		o.Xmin = utl.Min(o.Xmin, v.C[0])
		o.Xmax = utl.Max(o.Xmax, v.C[0])
		o.Ymin = utl.Min(o.Ymin, v.C[1])
		o.Ymax = utl.Max(o.Ymax, v.C[1])
		if nd > 2 {
			o.Zmin = utl.Min(o.Zmin, v.C[2])
			o.Zmax = utl.Max(o.Zmax, v.C[2])
		}

		// control points to initialise NURBS
		if has_nurbs {
			controlpts[i] = make([]float64, 4)
			for j := 0; j < 4; j++ {
				controlpts[i][j] = v.C[j]
			}
		}
	}

	// allocate NURBSs
	o.PtNurbs = make([]*gm.Nurbs, len(o.Nurbss))
	o.NrbFaces = make([][]*gm.Nurbs, len(o.Nurbss))
	for i, d := range o.Nurbss {
		o.PtNurbs[i] = new(gm.Nurbs)
		o.PtNurbs[i].Init(d.Gnd, d.Ords, d.Knots)
//.........这里部分代码省略.........

// testing_compare_results_u compares results with u-formulation
func TestingCompareResultsU(tst *testing.T, simfilepath, cmpfname, alias string, tolK, tolu, tols float64, skipK, verbose bool) {

	// FEM structure
	fem := NewFEM(simfilepath, alias, false, false, true, false, verbose, 0)

	// set stage
	err := fem.SetStage(0)
	if err != nil {
		chk.Panic("cannot set stage:\n%v", err)
	}

	// zero solution
	err = fem.ZeroStage(0, true)
	if err != nil {
		chk.Panic("cannot zero stage data:\n%v", err)
	}

	// read file with comparison results
	buf, err := io.ReadFile(cmpfname)
	if err != nil {
		tst.Errorf("TestingCompareResultsU: ReadFile failed\n")
		return
	}

	// unmarshal json
	var cmp_set T_results_set
	err = json.Unmarshal(buf, &cmp_set)
	if err != nil {
		tst.Errorf("TestingCompareResultsU: Unmarshal failed\n")
		return
	}

	// run comparisons
	dom := fem.Domains[0]
	dmult := 1.0
	for idx, cmp := range cmp_set {

		// displacements multiplier
		if idx == 0 && math.Abs(cmp.DispMult) > 1e-10 {
			dmult = cmp.DispMult
		}

		// time index
		tidx := idx + 1
		if verbose {
			io.PfYel("\n\ntidx = %d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n", tidx)
		}

		// load gofem results
		err = dom.Read(fem.Summary, tidx, 0, true)
		if err != nil {
			chk.Panic("cannot read 'gofem' results:\n%v", err)
		}
		if verbose {
			io.Pfyel("time = %v\n", dom.Sol.T)
		}

		// check K matrices
		if !skipK {
			if verbose {
				io.Pfgreen(". . . checking K matrices . . .\n")
			}
			for eid, Ksg := range cmp.Kmats {
				if e, ok := dom.Elems[eid].(*ElemU); ok {
					err = e.AddToKb(dom.Kb, dom.Sol, true)
					if err != nil {
						chk.Panic("TestingCompareResultsU: AddToKb failed\n")
					}
					chk.Matrix(tst, io.Sf("K%d", eid), tolK, e.K, Ksg)
				}
				if e, ok := dom.Elems[eid].(*Rod); ok {
					err = e.AddToKb(dom.Kb, dom.Sol, true)
					if err != nil {
						chk.Panic("TestingCompareResultsU: AddToKb failed\n")
					}
					chk.Matrix(tst, io.Sf("K%d", eid), tolK, e.K, Ksg)
				}
				if e, ok := dom.Elems[eid].(*ElastRod); ok {
					err = e.AddToKb(dom.Kb, dom.Sol, true)
					if err != nil {
						chk.Panic("TestingCompareResultsU: AddToKb failed\n")
					}
					chk.Matrix(tst, io.Sf("K%d", eid), tolK, e.K, Ksg)
				}
			}
		}

		// check displacements
		if verbose {
			io.Pfgreen(". . . checking displacements . . .\n")
		}
		for nid, usg := range cmp.Disp {
			ix := dom.Vid2node[nid].Dofs[0].Eq
			iy := dom.Vid2node[nid].Dofs[1].Eq
			chk.AnaNum(tst, "ux", tolu, dom.Sol.Y[ix], usg[0]*dmult, verbose)
			chk.AnaNum(tst, "uy", tolu, dom.Sol.Y[iy], usg[1]*dmult, verbose)
			if len(usg) == 3 {
				iz := dom.Vid2node[nid].Dofs[2].Eq
				chk.AnaNum(tst, "uz", tolu, dom.Sol.Y[iz], usg[2]*dmult, verbose)
//.........这里部分代码省略.........

func test_rwstep01(tst *testing.T) {

	verbose()
	chk.PrintTitle("rwstep01")

	buf, err := io.ReadFile("rw/data/beadpanel.step")
	if err != nil {
		tst.Errorf("cannot read file:\n%v", err)
		return
	}
	dat := string(buf)

	var stp rw.STEP
	err = stp.ParseDATA(dat)
	if err != nil {
		tst.Errorf("Parse filed:\n%v", err)
		return
	}

	var bsplines []*Bspline

	for _, scurve := range stp.Scurves {
		curve := stp.BScurves[scurve.Curve_3d]
		if curve == nil {
			continue
		}

		// collect vertices
		nv := len(curve.Control_points_list)
		verts := utl.DblsAlloc(nv, 4)
		for i, key := range curve.Control_points_list {
			if p, ok := stp.Points[key]; ok {
				for j := 0; j < 3; j++ {
					verts[i][j] = p.Coordinates[j]
				}
				verts[i][3] = 1.0
			} else {
				chk.Panic("cannot find point %q", key)
			}
		}

		// collect knots
		nk := 0
		for _, m := range curve.Knot_multiplicities {
			nk += m
		}
		knots := make([]float64, nk)
		k := 0
		for i, u := range curve.Knots {
			m := curve.Knot_multiplicities[i]
			for j := 0; j < m; j++ {
				knots[k] = u
				k++
			}
		}

		// create B-spline
		bsp := new(Bspline)
		bsp.Init(knots, curve.Degree)
		bsp.SetControl(verts)
		bsplines = append(bsplines, bsp)
	}

	if true {
		io.Pforan("n = %v\n", len(bsplines))
		for i, bsp := range bsplines {
			bsp.Draw3d("", "", 21, i == 0)
		}
		//plt.Show()
	}
}